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ABSTRACT 

This  study  investigates  and  develops  various  modifications  to  the  Multiple  Model  Adaptive 
Estimation  (MMAE)  algorithm.  The  standard  MMAE  uses  a  bank  of  Kalman  filters,  each  based  on 
a  different  model  of  the  system.  Each  of  the  filters  predict  the  system  response,  based  on  its  system 
model,  to  a  given  input  and  form  the  residual  difference  between  the  prediction  and  sensor 
measurements  of  the  system  response.  Model  differences  in  the  input  matrix,  output  matrix,  and 
state  transition  matrix,  which  respectively  correspond  to  an  actuator  failure,  sensor  failure,  and  an 
incorrectly  modeled  flight  condition  for  a  flight  control  failure  application,  were  investigated  in  this 
research.  An  alternative  filter  bank  structure  is  developed  that  uses  a  linear  transform  on  the  residual 
from  a  single  Kalman  filter  to  produce  the  equivalent  residuals  of  the  other  Kalman  filters  in  the 
standard  MMAE.  A  Neyman-Pearson  based  hypothesis  testing  algorithm  is  developed  that  results  in 
significant  improvement  in  failure  detection  performance  when  compared  to  the  standard  hypothesis 
testing  algorithm.  Hypothesis  testing  using  spectral  estimation  techniques  is  also  developed  which 
provides  superior  failure  identification  performance  at  extremely  small  input  levels. 
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PRACTICAL  IMPLEMENTATION  OF 


MULTIPLE  MODEL  ADAPTIVE  ESTIMATION  USING 
NEYMAN-PEARSON  BASED  HYPOTHESIS  TESTING  AND 
SPECTRAL  ESTIMATION  TOOLS 


I.  Introduction 


1.1  Chapter  Overview 

In  this  chapter  we  introduce  the  development  of  modifications  to  the  Multiple  Model 
Adaptive  Estimation  (MMAE)  algorithm.  To  demonstrate  the  MMAE  performance  with  these 
modifications,  we  have  chosen  to  use  a  flight  control  failure  detection  application,  which  is  described 
in  Section  1.2.  We  also  will  point  out  some  of  the  difficulties  encountered  with  contemporary 
MMAE  designs  that  motivated  the  development  of  these  modifications.  Finally,  in  Section  1.3  we 
present  an  overview  of  the  structure  of  this  dissertation. 

1.2  Description  of  the  Motivating  Problem 

The  MMAE  algorithm  shows  tremendous  promise  in  enhancing  flight  control  performance. 
The  constant  pressure  to  improve  aircraft  performance,  particularly  in  military  aircraft,  has  forced 
designers  to  build  aircraft  with  instabilities  that  are  beyond  that  capabilities  of  human  pilots  to 
control.  Thus,  these  designers  have  been  forced  to  rely  on  increasingly  complex  flight  control 
systems.  These  systems  are  designed  for  a  particular  system  configuration  (i.e.  a  fully  functional 
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Figure  1.  Multiple  Model  Adaptive  Estimation  Algorithm 


airaaft)  and  usually  cannot  perform  adequately  when  the  system  is  configured  substantially 
differently  (i.e.  a  flight  control  failure).  Designers  are  constructing  several  flight  control  system 
designs  that  will  perform  adequately  for  different  system  configurations,  but  the  flight  control  system 
must  be  informed  which  configuration  is  the  appropriate  one  to  choose.  One  such  method  of 
detecting  and  identifying  system  failures  is  the  MMAE. 

The  MMAE,  diagrammatically  shown  in  Figure  1,  is  composed  of  a  bank  of  Kalman  filters, 
each  imbued  with  its  own  model,  and  a  hypothesis  testing  algorithm.  Each  elemental  Kalman  filter 
uses  its  own  model,  along  with  a  given  input  («),  to  develop  an  estimation  of  the  current  aircraft 
states  (jCjfc),  independent  of  the  other  filters.  The  filter  then  uses  this  estimate,  along  with  the  current 
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measurement  of  those  states  (z),  to  form  the  residual  (r^),  which  is  the  difference  between  the 
measurement  and  the  filter's  prediction  of  the  measurements  before  they  arrive.  The  residuals  from 
the  filters  are  used  by  the  hypothesis  testing  algorithm  as  a  relative  indication  of  how  close  each  of 
the  filter  models  are  to  the  true  model,  i.e.,  to  the  real-world  situation.  The  smaller  the  residual,  the 
closer  the  filter  model  matches  the  true  model.  The  hypothesis  testing  algorithm  first  scales  the 
residuals  to  account  for  various  uncertainties  and  noises  in  the  measurements  (as  developed  in  detail 
in  Section  3.3.1),  and  then  computes  the  conditional  probability  for  each  of  the  hypotheses  modeled 
in  the  bank  of  Kalman  filters  (p*).  These  probabilities  are  then  used  to  weight  the  individual  Kalman 
filter  state  estimates  to  produce  a  blended  estimate  of  the  true  states  (Xmmae)’  which  can  then  be  used 
as  the  optimal  estimate  of  the  states  by  a  control  sj^tem.  When  used  for  failure  identification,  each 
of  the  Kalman  filters  would  model  a  different  failure  condition  and  the  residuals  from  each  filter 
would  indicate  how  close  that  filter's  model  is  to  the  actual  failure  condition.  By  monitoring  these 
residuals,  the  hypothesis  tester  can  estimate  the  current  failure  status  of  the  airaaft  (a). 

The  specific  aircraft  model  [62]  that  is  used  for  this  application  was  developed  for  the 
LAMBDA  flight  vehicle.  The  LAMBDA  is  an  unmanned  research  vehicle  developed  by  the  Flight 
Control  Division  of  the  Flight  Dynamics  Directorate,  Wright  Laboratory,  as  an  affordable,  flexible 
research  vehicle  for  testing  and  demonstrating  flight  control  concepts,  devices,  and  systems  [62]. 

The  flight  condition  of  the  LAMBDA  is  determined  by  five  parameters:  the  aircraft  weight,  forward 
velocity,  dynamic  pressure,  center  of  gravity,  and  trim  angle.  The  average  flight  condition 
parameters  are  a  weight  of  200  lbs,  a  speed  of  160  ft/sec,  a  dynamic  pressure  of  30.43  Ib/fF,  a  center 
of  gravity  located  at  46.8  inches  from  the  aircraft  nose,  and  a  trim  angle  of  zero  degrees.  The 
LAMBDA  usually  flies  at  low  altitude;  therefore  the  dynamic  pressure  ('/2pv^)  is  directly  related  to 
forward  velocity  (v)  since  the  air  density  (p)  is  essentially  constant.  The  trim  angle  is  usually  very 
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Figure  2.  LAMBDA  Flight  Envelope 


small  for  a  well  trimmed  airaaft;  therefore  we  assume  a  trim  angle  of  zero.  Since  the  altitude  and 
trim  angle  are  considered  constant,  we  end  up  with  a  three-dimensional  flight  envelope  for  the 
LAMBDA  as  shown  in  Figure  2,  which  also  shows  the  design  flight  condition.  This  flight  condition 
is  considered  average  because  the  LAMBDA  is  usually  flown  close  to  maximum  velocity,  at  low 
altitude,  with  the  variation  in  weight  primarily  dependent  on  the  amount  of  fuel  consumed.  Therefore 
the  design  flight  condition  is  a  median  value  of  the  normal  operating  points. 

Embedded  within  the  linearized  model  (linear  perturbation  model  about  trim  condition)  of 
the  LAMBDA  are  several  other  assumptions  that  were  used  to  develop  the  model  [62],  namely  lateral 
and  longitudinal  axis  decoupling  and  first  order  actuator  dynamics.  The  linearized  model  was 
developed  by  estimating  flight  coefficients  from  geomeu-ic  data  and  then  adjusting  these  estimates 
using  good  engineering  judgment  and  flight  test  data.  For  a  fully  functional  aircraft,  the  cross-axis 
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coupling  terms  are  so  small  that  they  produce  a  negligible  airaaft  response.  Currently,  all  test  flights 
have  been  conducted  using  a  fully  functional  aircraft  [62],  and  therefore,  the  flight  test  data  does  not 
contain  enough  information  to  estimate  the  aoss-axis  terms.  The  assumption  of  first  order  actuator 
dynamics  compares  quite  well  with  the  flight  test  data. 

To  model  single  actuator  failures,  we  have  assumed  that  a  single  flight  control  surface  failure 
will  produce  half  of  the  expected  response  from  dual  control  surfaces.  For  instance,  a  certain 
elevator  input  might  produce  a  10  degree  pitch-up  response.  We  have  assumed  that,  for  this  same 
input,  a  failure  of  only  the  right  elevator  actuator  would  produce  half  of  this  response  (5  degree  pitch- 
up)  since  only  the  left  elevator  actuator  would  be  functional,  thereby  decreasing  the  actuator  surface 
area  by  one  half. 

As  mentioned  above,  we  have  assumed  that  the  cross-axis  coupling  terms  are  negligible  even 
for  single  surface  failures,  primarily  because  there  is  insufficient  data  to  estimate  these  terms  [62]. 
Using  our  example,  a  right  elevator  failure  would  produce  a  small  yawing  and  rolling  response, 
which  we  assume  to  be  negligible.  An  exception  to  this  assumption  is  the  cross-coupling  between 
the  ailerons  and  the  yaw  axis  because  the  long  yaw  moment  arm  of  the  ailerons  could  easily  produce 
a  significant  yaw.  Fortunately,  there  is  sufficient  data  to  estimate  the  aileron's  effect  on  the  yaw  rate 
(r)  and  the  rudder's  effect  on  the  roll  rate  (p).  At  the  design  point,  the  contribution  of  the  aileron  to  a 
yaw  rate  is  10.5%  and  the  contribution  of  the  rudder  to  a  roll  rate  is  only  0.3%.  This  supports  our 
assumption  that  the  cross-axis  terms  are  small.  However,  Swift  found  noticeable  roll/pitch  coupling 
dynamics  for  the  LAMBDA,  but  was  unable  to  estimate  the  cross  coupling  terms  [62:  3.50].  A  flight 
test  with  single  actuator  failures  should  provide  the  much  needed  data  for  estimating  these  cross-axis 
coupling  terms,  and  the  aircraft  model  can  then  be  corrected  for  further  tests. 
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Previous  research  [19, 31,41]  that  developed  an  MMAE-based  failure  detection  algorithm 
for  this  flight  vehicle  identified  three  problems  with  contemporary  MMAE  designs.  These  problems 
are  a  lack  of  a  methodical  way  to  make  a  definitive  failure  status  declaration,  the  computational  cost 
for  implementing  these  designs,  and  the  need  to  dither  the  system  constantly  to  provide  adequate 
excitation  for  good  failure  detection.  These  problems  are  described  in  more  detail  below. 

Most  flight  control  systems  require  a  definite  failure  status  declaration  from  the  failure 
detection  and  identification  algorithm.  The  conditional  probabilities  produced  by  the  MMAE  give 
only  a  relative  indication  of  the  failure  status.  To  make  a  declaration  of  the  failure  status, 
contemporary  MMAE  designers  choose  a  probability  threshold  that  is  used  to  declare  if  a  certain 
hypothesized  failure  has  occurred.  If  the  conditional  probability  for  that  hypothesis  exceeds  the 
threshold,  then  the  hypothesized  failure  would  be  declared.  Unfortunately,  there  is  no  method  of 
choosing  this  threshold  aside  from  extensive  computer  simulations.  These  simulations  would  give  a 
relative  indication  of  the  tradeoff  between  the  false  alarm  rate  and  failure  detection  performance  and 
the  designer  would  use  good  engineering  judgement  to  choose  the  threshold.  A  more  satisfactory 
design  method  would  be  to  be  able  to  dictate  the  maximum  acceptable  false  alarm  rate  and  desired 
failure  detection  probability,  and  have  these  two  parameters  dictate  the  threshold.  Such  a  technique 
is  available  using  a  Neyman-Pearson  base  hypothesis  testing  algorithm  that  will  be  developed  in  this 
research. 

MMAE  can  be  a  very  costly  algorithm  to  implement.  The  Kalman  filter  algorithm, 
presented  in  Section  3.2.1,  can  be  computationally  intensive,  particularly  if  a  large  number  of  filter 
states  is  needed  to  capture  the  real-world  aircraft  dynamics.  Using  a  bank  of  Kalman  filters  further 
compounds  the  computational  cost  of  implementing  an  MMAE  algorithm.  MMAE  can  be 
implemented  on  a  set  of  parallel  processors,  because  of  its  inherent  parallel  structure,  but  the 
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required  power,  cooling,  and  data  bus  interconnection,  along  with  the  weight  and  cost,  of  producing 
this  set  of  parallel  processors  could  easily  make  this  implementation  infeasible.  We  develop  an 
equivalent  MMAE  filter  bank  structure  that  uses  the  residual  and  state  estimates  from  a  single 
Kalman  filter,  along  with  the  linear  transformation  that  is  developed  in  Section  3.2.5,  to  produce  the 
equivalent  residual  from  another  Kalman  filter.  This  transformation  computes  the  equivalent 
residual  by  using  the  known  differences  between  the  two  Kalman  filter  models.  These  model 
differences  usually  produce  sparse  matrices,  which  significantly  reduces  the  computational  cost  of 
computing  the  equivalent  residuals.  Thus,  the  bank  of  filters  can  be  replaced  by  a  single  filter  with 
several  linear  transformations,  each  of  which  produces  the  equivalent  residual  to  the  Kalman  filter 
that  it  replaces,  but  with  a  reduction  in  the  cost  of  computing  the  residual. 

The  primary  objection  to  implementing  an  MMAE-based  failure  detection  algorithm  is  the 
need  to  dither  the  system  constantly.  The  MMAE  compares  the  magnitudes  of  the  residuals 
(appropriately  scaled  to  account  for  various  uncertainties  and  noises)  from  the  various  filters  and 
chooses  the  hypothesis  that  corresponds  to  the  residual  that  has  a  history  of  having  smallest  (scaled) 
magnitude.  Large  residuals  need  to  be  produced  by  the  filters  with  models  that  are  incorrect  to  be 
able  to  identify  these  incorrect  hypotheses.  The  residual  is  the  difference  between  the  measurement 
of  the  system  output  and  the  filter's  prediction  of  what  that  measurement  should  be,  based  on  the 
filter-assumed  system  model.  Therefore,  to  produce  the  needed  large  residuals  in  the  incorrect  filters, 
we  need  to  produce  a  history  of  sufficiently  large  system  outputs,  so  we  need  to  dither  the  system 
constantly  and  thereby  excite  the  system  states.  For  flight  control  failure  applications,  we  would 
need  to  move  the  airaaft  continually  in  all  axes  to  produce  the  desired  failure  detection  performance, 
to  which  most  pilots  (not  to  mention  passengers)  would  strenuously  object.  We  develop  a  modified 
MMAE  algorithm  that  uses  the  correlation  of  the  residual  with  signal  processing  techniques  to 
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produce  the  desired  failure  detection,  while  significantly  reducing  the  objectionable  dither,  possibly 
to  the  point  of  being  subliminal. 

1.3  Dissertation  Overview 

In  the  subsequent  chapters  we  develop  this  research  in  more  detail.  We  have  presented  a 
general  introduction  to  the  MMAE  and  a  brief  description  of  the  problems  that  have  motivated  this 
research.  Chapter  2  gives  a  brief  background  on  the  development  of  the  MMAE  and  a  description  of 
some  statistical  signal  processing  techniques  that  are  used  in  the  research.  Chapter  3  presents  the 
development  of  the  theory  behind  the  various  modifications  to  the  MMAE  that  have  been  developed. 
Chapter  4  presents  the  results  from  a  specific  simulation  of  these  modifications  using  the  flight 
control  application  described  above.  Finally,  in  Chapter  5  we  present  a  comparison  of  the  various 
modifications,  draw  pertinent  conclusions,  and  finish  with  our  recommendations  for  future  research. 
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n.  Background 


2. 1  Chapter  Overview 

In  this  chapter  we  introduce  the  background  research  that  has  led  up  to  this  investigation. 

We  start,  in  Section  2.2,  with  a  brief  survey  of  the  development  of  the  Multiple  Model  Adaptive 
Estimation  (MMAE)  and  describe  some  of  the  research  of  its  use  for  various  applications.  In  Section 
2.3  we  present  an  overview  of  statistical  signal  processing  techniques  that  may  help  enhance  the 
performance  of  MMAE.  These  techniques  include  hypothesis  testing  based  on  various  spectral 
estimation  techniques  and  the  Neyman-Pearson  lemma . 

2.2  Multiple  Model  Adaptive  Estimation 

2.2.1  Multiple  Model  Adaptive  Estimation  Development.  The  use  of 

multiple  filters  in  a  parallel  structure  to  generate  adaptive  estimation  algorithms  was  first  developed 
by  D.  T.  Magill  [35].  He  arranged  a  number  of  Kalman  filters,  each  with  different  time  invariant 
plant  models,  and  used  the  residuals  from  these  filters  to  form  an  appropriately  weighted  sum  of  the 
Kalman  filter  estimates,  as  shown  in  Figure  1 .  He  showed  that  this  adaptive  estimation  algorithm 
produced  the  optimal  estimate  in  the  minimum  mean  square  error  sense  for  a  Gauss-Markov  process, 
over  a  finite  set  of  models,  one  of  which  is  the  "true"  model  which  accurately  represents  the  true 
system.  This  work  was  extended  to  handle  equivalent  discrete-time  systems  by  C.  B.  Chang  and  M. 
Athans  [6,  7,  8]. 

The  properties  of  this  multiple  model  adaptive  estimation  concept  were  further  developed  by 
Lainiotis  and  others  [10,  20,  28,  29].  They  investigated  the  properties  of  this  concept  when  using  a 
discrete  number  of  models  to  represent  a  continuous  domain  of  plant  models  [28, 29].  With  others. 
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he  showed  that  good  performance  was  attained  with  tightly  tuned  Kalman  filters,  but  these  filters 
produced  erratic  behavior  if  the  filter  with  the  correct  model  was  not  in  the  filter  bank  [10, 20]. 

A.  S.  Willsky  [68]  surveyed  a  number  of  failure  detection  methods,  using  a  performance 
index  as  an  empirical  measurement  of  the  capabilities  of  the  various  methods.  He  surveyed  three 
classes  of  failure  detection  methods.  The  first  was  the  use  of  specific  failure-sensitive  filters,  such  as 
limited  memory  or  Fagin  age  weighting  filter  and  filters  that  included  failure  states  in  the  filter 
model.  A  second  class  was  methods  that  conduct  statistical  tests  on  the  filter  innovations,  such  as  the 
Generalized  Likelihood  Ratio  (GLR)  and  voting  techniques.  The  third  class  consisted  of  multiple 
hypothesis  filter  detectors,  specifically  the  MMAE.  The  performance  index  he  used  included  the 
types  of  failure  modes  that  can  be  detected,  implementation  complexity,  various  performance 
measures  (false  alarms,  detection  delays,  repeatability),  and  robustness  in  the  presence  of  modeling 
errors.  He  found  that  the  multiple  model  adaptive  estimator  will  yield  the  best  performance  over  the 
widest  class  of  failures. 

2.2.2  Multiple  Model  Adaptive  Control  Development.  M.  Athans  et.  al., 

extended  this  adaptive  estimation  concept  to  adaptive  control  developed  for  the  flight  control  system 
of  NASA's  F-8C  flight  test  airaaft  [1].  They  weighted  the  optimal  control  signals  generated  by  a 
bank  of  Linear  system-Quadratic  cost  function-Gaussian  noise  distribution  (LQG)  controllers,  each 
with  an  embedded  Kalman  filter  using  a  different  plant  model,  and  then  summed  these  weighted 
control  signals  together,  as  shown  in  Figure  3.  They  found  that  this  algorithm  provided  good  control 
at  flight  conditions  that  were  close  to  the  design  conditions  of  the  Kalman  filter  models,  as  long  as 
the  system  was  appropriately  excited,  using  a  test  signal  or  dither,  to  attain  good  failure 
identification.  They  also  found  that  this  algorithm  was  sensitive  to  high  frequency  noise,  specifically 
strong  wind  gusts,  and  required  low  pass  filtering  to  attain  good  performance. 
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Figure  3.  Multiple  Model  Adaptive  Control  Algorithm 


Greene  and  Willsky  [16]  examined  and  defined  stability  regions  where  the  multiple  model 
adaptive  control  (MMAC)  algorithm  yields  non-oscillatory  responses.  They  found  that  the  stability 
of  the  MMAC  was  determined  by  the  relation  between  the  growth  rate  of  the  most  imstable  mode  of 
a  Kalman  filter  with  a  mismatched  model,  when  compared  to  the  truth  model,  and  the  rate  of  decay  of 
the  slowest  stable  mode  of  a  Kalman  filter  with  a  matching  model.  They  presented  a  method  of 
computing  the  borders  where  the  algorithm  is  neutrally  stable,  which  then  defines  the  "domain  of 
attraction,"  or  the  region  where  the  MM  AC's  response  does  not  oscillate. 

Longitudinal  control  of  the  Short  Take-Off  and  Landing  (STOL)  F-15  using  Multiple  Model 
Adaptive  Control  (MMAC)  was  investigated  by  Pogoda  [42, 52]  and  Stevens  [43,  60].  Pogoda 
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developed  an  algorithm  for  the  landing  phase  of  the  flight  regime,  which  would  reconfigure  the  flight 
control  in  the  presence  of  a  single  failed  control  surface  or  sensor.  He  designed  Kalman  filter  models 
and  controllers  for  a  fully  functional  aircraft,  a  failed  stabilator,  a  failed  "pseudo-surface,"  and  a 
failed  pitch  rate  sensor.  The  "pseudo-surface"  was  a  combination  of  the  canards,  ailerons,  and  flaps, 
to  allow  a  reduction  in  the  number  of  aircraft  states  and  independent  control  surfaces.  Stevens 
extended  this  research  further  by  including  failed  reverser  vanes,  a  failed  velocity  sensor,  and  a  failed 
flight  path  sensor.  He  first  investigated  "soft"  or  partial  failures  of  these  sensors  and  surfaces,  with 
the  soft  failures  modeled  as  either  partial  power  to  a  flight  surface,  an  increase  in  the  sensor  noise,  or 
an  increase  in  the  sensor  bias.  He  then  investigated  the  performance  of  a  hierarchical  structure  of 
multiple  model  adaptive  controllers  to  detect  the  presence  of  double  failures,  such  as  a  stabilator 
surface  failure  followed  by  a  flight  path  angle  sensor  failure.  This  hierarchical  structure  of 
controllers  started  with  a  bank  of  controllers,  with  different  single  failure  models  (including,  of 
course,  a  no-failure  model),  that  would  detect  the  first  failure  and  then  switch  to  another  bank  of 
controllers  that  had  models  that  assumed  both  the  detected  failure  and  a  second  failure  (including  no 
second  failure).  Each  of  the  banks  also  had  a  controller  with  a  failure  model  that  assumed  that  the 
detected  failure  had  not  occurred  after  all,  which  allowed  the  structure  to  correct  any 
misidentifications.  Pogoda  and  Stevens  both  found  that  the  MM  AC  structure  was  able  to  identify  the 
failures  properly  and  reconfigure  the  control  law  to  maintain  stable  flight  control.  Stevens  found  that 
the  MMAC  would  blend  the  appropriate  controller  commands  in  the  presence  of  soft  failures,  and 
that  the  hierarchical  structure  would  properly  detect  multiple  failures. 

2,2.3  Multiple  Model  Adaptive  Estimation  Based  Control  in  a  similar 

manner,  Stratton  [61]  and  Menke  [47, 48, 49]  have  developed  an  MMAE-based  control  algorithm 
for  both  the  longitudinal  and  lateral  axes  of  the  VISTA  F-16  airaaft.  In  contrast  to  the  MMAC  of 
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Figure  3,  this  algorithm  uses  the  MMAE  of  Figure  1,  feeding  and  ®  single  controller  block; 

in  their  implementation,  only  £,aiAK  ^  ^  ^  already  designed  and  validated  flight 

control  system.  Stratton  used  a  flight  condition  of  Mach  0.8  and  altitude  10,000  feet,  while  Menke 
designed  for  a  flight  condition  of  Mach  0.4  and  altitude  20,000  feet,  the  latter  involving  low  dynamic 
pressures  and  thus  presenting  a  more  difficult  failure  detection  problem.  They  investigated  both 
single  and  multiple,  hard  and  soft,  actuator  and  sensor  failures.  Included  in  their  study  was  the  effect 
of  a  test  signal  or  commanded  dither  to  aid  in  identifying  the  failure  during  benign  straight  and  level 
flight  conditions.  Several  dither  signals  were  tested,  including  sine  waves,  square  waves,  triangular 
waves,  and  pulse  trains,  at  levels  that  were  deemed  either  subliminal  (up  to  ±  0.1  g's  in  the 
longitudinal  axis  and  ±  0.2  g's  in  the  lateral  axis)  or  non-subliminal  (reasonable  physical  acceleration 
limits  at  the  pilot's  station).  They  found  that  the  MMAE  identified  these  failures  correctly  as  long  as 
an  appropriate  dither  signal  was  present  that  excited  all  failure  modes  in  both  axes. 

Hide  [11, 12]  and  Stepaniak  [58, 59]  continued  the  development  of  the  MMAE-based 
control  algorithm  for  the  VISTA  F-16.  Hide  tested  the  algorithm  against  a  full-scale  nonlinear  truth 
model,  instead  of  a  linearized  truth  model  that  was  used  for  the  previous  performance  analyses. 
Initially,  he  found  fairly  poor  failure  identification  performance  that  was  attributed  to  a  mismatch 
between  the  full-scale  nonlinear  model  and  the  linearized  design  model.  He  tuned  the  Kalman  filters 
and  achieved  much  better  performance,  thus  showing  the  robustness  of  the  MMAE  algorithm,  even 
when  operating  with  a  nonlinear  system  while  the  MMAE  Kalman  filters  are  based  on  linearized 
models.  Stepaniak  used  an  MMAE  to  redistribute  the  control  input  to  redundant  nonfailed  actuators 
using  the  existing  VISTA  F-16  control  system.  He  obtained  excellent  results,  showing  that  this 
MMAE-based  control  algorithm  could  track  a  desired  state  trajectory,  as  long  as  the  redundant 
actuators  were  not  forced  into  saturation.  Thus,  he  showed  that  the  MMAE-based  control  algorithm 
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can  completely  compensate  for  a  failed  actuator,  up  to  the  limits  of  command  authority,  if  redundant 
control  surfaces  are  available. 

2.2.4  Multiple  Model  Adaptation  Applications.  Multiple  model  adaptive 

algorithms  have  been  successfully  developed  for  a  number  of  other  applications.  It  was  investigated 
for  the  detection  and  tracking  of  maneuvering  targets  [8, 14, 15,  37, 44,  50, 51,  63, 64,  65],  flexible 
space  structure  control  [13, 23, 24, 32, 33, 54],  multiple  hypotheses  testing  [1, 45],  and  prevention 
of  the  initial  divergence  of  extended  Kalman  filters  due  to  large  initial  uncertainties  [22, 46].  It  has 
been  studied  for  use  in  diverse  applications  such  as  instrument  failure  detection  in  a  pressurized 
water  fusion  reactor  [10],  autonomous  monitoring  of  cardiac  patients  [17],  adaptive  signal 
processing  of  seismic  data  [58],  and  detection  of  incidents  on  freeways  [69]. 

2.3  Signal  Proce.s.sing  Techniques 

2.3.1  Spectral  Estimation  Techniques.  Kay  [26]  presents  a  methodical  survey  of 

some  of  the  many  spectral  estimation  techniques.  He  groups  these  techniques  into  classical  methods, 
those  that  are  based  on  Fourier  transform  and  filtering  theory,  and  modem  techniques  that  are  based 
on  time  series  analysis  and  filtering  theory.  In  the  classical  group  he  describes  periodograms, 
averaged  periodograms,  and  Blackman-Tukey  spectral  estimators.  Based  on  the  definition  of  the 
power  spectral  density  of  a  signal,  the  Fourier  transform  of  the  autocorrelation  function,  these 
techniques  window  the  data  (using  various  shaped  windows  for  each  of  the  different  techniques)  and 
then  compute  the  spectral  content  across  the  data  window  at  various  frequencies  by  computing  the 
Fourier  transform  of  an  estimate  of  the  autocorrelation  function  of  the  windowed  data.  Modem 
techniques  are  based  on  parametric  modeling  of  the  signal,  with  the  assumption  that  the  signal  is 
composed  of  a  known  number  of  sinusoids. 
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The  primary  motivation  for  the  research  of  these  techniques  is  detecting  and  characterizing 
sinusoids  with  time-varying  amplitude,  embedded  in  a  background  of  noise,  which  describes  many 
applications  where  the  data  set  is  relatively  short  compared  to  the  duration  of  the  sinusoid.  This 
occurs  with  speech  synthesis,  where  a  speech  sound  may  only  last  for  about  20  to  80  msec  [36],  and 
in  Doppler  radar  and  sonar,  where  the  propagation  characteristics  of  the  transition  medium  may 
change  in  time  [25] .  In  other  applications  the  small  data  set  is  due  to  a  genuine  lack  of  data. 

Seismic  data  sets  are  transient  because  events  such  as  volcanic  eruptions  and  earthquakes  last  for 
very  short  periods  of  time  [30].  Other  applications  have  small  data  sets  because  of  the  prohibitive 
cost  of  collecting  the  data,  such  as  optical  interferometry  [4],  vibration  analysis  [2],  radio  astronomy 
[67],  image  processing  [18],  and  the  application  that  we  are  studying,  flight  control  failure  detection. 
There  are  many  other  methods  of  signal  analysis  that  apply,  such  as  Maximum  Entropy  Methods 
(MEM),  particularly  useful  for  small  data  sets  [53],  but  we  have  chosen  to  research  the  use  of 
periodograms  (presented  in  Section  3.2.6)  because  it  has  been  extensively  researched,  is 
straightforward  to  implement,  and  provides  a  basic  conceptual  framework  for  implementing  the  other 
techniques. 

2.3.2  Nevman-Pear.son  Hypothesis  Te.sting.  Scharf  [53:103-166  ]  presents  the 

development  of  a  Neyman-Pearson  Detector.  The  basis  of  the  detector  is  a  hypothesis  test  in  which 
the  distribution  of  a  test  statistic  is  known,  for  each  of  the  hypotheses  that  is  being  tested.  The 
hypothesis  test  is  based  on  the  Neyman-Pearson  lemma,  which  will  be  formally  presented  in  Section 
3.3.2,  which  shows  that  a  likelihood  ratio  that  is  formed  using  these  known  distributions  of  the  test 
statistic,  provides  the  best  probability  of  detection  for  a  given  false  alarm  probability,  for  a  fixed 
number  of  data  samples.  The  probability  of  detection  is  defined  as  the  probability  that  the  correct 
hypothesis  is  actually  chosen,  given  that  the  incorrect  hypothesis  is  considered  to  be  in  force  at  the 
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time  of  the  hypothesis  test.  For  example,  the  probability  of  detection  is  the  probability  that,  if  a 
failure  occurs,  the  hypothesis  test  chooses  the  failure  hypothesis  when  prior  to  the  test  the  no  failure 
hypothesis  was  correct.  The  false  alarm  probability  is  defined  as  the  probability  that  the  incorrect 
hypothesis  is  chosen,  given  that  the  correct  hypothesis  is  considered  in  force  at  the  time  of  the 
hypothesis  test.  For  example,  this  is  the  probability  that,  if  a  failure  does  not  occur,  the  hypothesis 
test  chooses  the  failure  hypothesis,  when  prior  to  the  test  the  no  failure  hypothesis  was  correct.  Thus, 
we  want  to  maximize  the  probability  of  detection,  while  minimizing  the  probability  of  false  alarm, 
which  illustrates  the  importance  of  the  Neyman-Pearson  lemma. 

These  definitions  will  differ  slightly  from  the  usual  convention  when  the  hypothesis  that  is 
considered  in  force  is  not  the  fully  functional  hypothesis.  To  illustrate  this,  assume  that  the  left 
elevator  has  failed  and  that  hypothesis  is  in  force.  If  a  hypothesis  test  is  conducted,  and  the  test 
incorrectly  chooses  the  fully  functional  hypothesis,  that  would  be  a  false  alarm  by  the  definitions 
given  above.  The  usual  convention  is  to  call  this  incorrect  declaration  a  missed  detection  because  the 
test  mistakenly  missed  the  detection  of  the  failed  elevator.  We  will  use  the  definitions  given  above  to 
be  consistent  with  the  development  of  the  Neyman-Pearson  hypothesis  test. 

Frequently,  the  fixed  length  of  the  data  set  will  not  provide  enough  information  to  be  able  to 
distinguish  the  various  hypotheses  from  each  other.  If  the  number  of  data  samples  is  not  set,  thus 
defining  a  sequential  estimation  problem,  the  sequential  probability  ratio  test  (SPRT)  is  superior  to 
the  Neyman-Pearson  most  powerful  test  because  it  will  not  make  a  decision  until  sufficient 
information  is  obtained  to  distinguish  the  hypotheses  from  each  other  [34].  We  have  chosen  to  use 
the  Neyman-Pearson  most  powerful  test  to  lay  the  conceptual  framework  for  this  type  of  hypothesis 
testing.  Scharf  [53]  also  develops  the  methodology  for  selecting  the  decision  threshold  that  is  used  to 
choose  between  the  various  hypotheses.  The  special  case  of  a  test  statistic  that  is  normally 


16 


distributed  with  differing  means  is  particularly  important  for  our  application,  since  the  Kalman  filter 
residual  will  be  shown,  in  Section  3.2,  to  be  normally  distributed  with  a  nonzero  mean  if  the  Kalman 
filter  model  is  inaccurate,  and  thus  we  can  use  the  residual  as  a  test  statistic. 


2.4  Chapter  Summary 

In  this  chapter  we  have  presented  a  brief  summary  of  the  research  that  has  led  up  to  this 
investigation.  First,  we  described  the  development  of  the  MMAE,  along  with  two  methods  of  using 
the  multiple  model  concept  for  control  applications.  The  MMAC  algorithm  uses  a  bank  of  the 
Kalman  filters  with  LQG  controllers  that  were  each  developed  for  a  specific  hypothesis,  and  then 
uses  the  conditional  probabilities  from  the  hypothesis  testing  algorithm,  to  weight  the  commanded 
inputs  from  the  various  controllers.  The  MMAE-based  control  algorithm  uses  an  MMAE  to  provide 

estimates,  rather  than  raw  measurements,  to  a  single  separate  control  system.  The  MMAE-based 

/ 

control  algorithm  with  control  redistribution  has  been  shown  to  be  able  to  compensate  completely  for 
any  failed  airaaft  sensor  and  for  any  failed  actuator,  up  to  the  point  of  saturation  of  the  actuator,  if 
redundant  control  surfaces  are  available.  We  briefly  listed  some  of  the  applications  that  have  used 
MMAE.  Then,  we  described  some  spectral  estimation  techniques  which  show  great  promise  in 
enhancing  the  failure  detection  performance  of  the  MMAE.  We  have  chosen  to  focus  on  the 
periodogram  spectral  estimation  technique  because  its  characteristics  are  well  researched  and  it 
provides  an  excellent  conceptual  framework  for  developing  other  spectral  estimation  techniques.  We 
introduced  the  Neyman-Pearson  detector,  which  we  will  use  as  a  framework  for  developing  an 
alternative  hypothesis  testing  algorithm. 


17 


HI.  Theory  Development 


3.1  Chapter  Overview 

A  Multiple  Model  Adaptive  Estimator  (MMAE)  consists  of  a  bank  of  parallel  Kalman 
filters,  each  with  a  different  model,  and  a  hypothesis  testing  algorithm  as  shown  in  Figure  4. 


Figure  4.  Multiple  Model  Adaptive  Estimation  Algorithm 


The  internal  models  of  the  Kalman  filters  can  be  represented  by  a  discrete  value  of  a  parameter  vector 
(fl*).  The  Kalman  filters  are  provided  a  measurement  vector  (z)  and  the  input  vector  («),  and 
produce  a  state  estimate  and  a  residual  (r*).  The  hypothesis  testing  algorithm  uses  the  residuals 
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to  compute  conditional  probabilities  (pi)  of  the  various  hypotheses  that  are  modeled  in  the  Kalman 
filters  conditioned  on  the  history  of  measurements  received  up  to  that  time,  and  an  estimate  of  the 
true  parameter  vector  (a).  For  the  flight  control  failure  application  in  which  we  are  interested,  each 
Kalman  filter  has  a  different  failure  model  (a,)  that  it  uses  to  form  the  state  estimate  (x^)  and  the 
residual  (r^).  The  hypothesis  testing  algorithm  produces  an  estimate  of  the  failure  status  of  the  flight 
control  system  (a).  The  Kalman  filter  equations  will  be  presented  in  Section  3.2.  In  Sections  3.2.2 
through  3.2.4,  we  will  show  that  the  magnitude  of  the  residual  from  a  particular  filter  can  be  used  as 
a  relative  measure  of  the  inaccuracy  (conversely,  the  accuracy)  of  the  failure  model  used  by  that  filter. 
In  Section  3.3  we  will  study  various  hypothesis  testing  algorithms  that  assign  conditional 
probabilities  (pi)  to  each  of  the  hypotheses  that  were  used  to  form  the  Kalman  filter  models.  These 
conditional  probabilities  indicate  the  relative  correcmess  of  the  various  filter  models,  and  can  be  used 
to  select  the  best  estimate  of  the  true  system  model,  weight  the  individual  state  estimates 
appropriately,  and  form  a  probability-weighted  average  state  estimate  (Jcmmae)-  T'he  specific 
workings  of  the  algorithms  in  these  blcKks,  along  with  various  modifications  to  those  algorithms,  are 
examined  in  the  following  sections. 

3.2  Kalman  Filter  Bank 

3.2.1  Basic  Equations.  The  Kalman  filter  algorithm  uses  a  model  of  the  true  system 

to  generate  estimates  of  the  noise-corrupted  measurements  before  they  are  taken  at  a  particular 
sample  time,  assuming  that  the  measurements  are  linear  combinations  of  the  system  states.  Previous 
researchers  have  computed  the  mean  and  covariance  of  the  residual  under  the  assumption  that  the 
Kalman  filter  model  accurately  represents  the  true  system.  We  will  relax  this  assumption  by  allowing 
certain  differences  between  the  Kalman  filter  model  and  a  model  that  accurately  represents  the  true 
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system.  Under  these  different  assumptions,  we  will  compute  the  mean  and  covariance  of  the  Kalman 
filter  residual. 

We  assume  that  we  have  noise-corrupted  measurements  of  a  possibly  nonlinear  system  that 
we  are  examining.  These  true  system  measurements  will  be  denoted  z(t,).  We  also  assume  that  there 
exists  a  linear  model  that  will  produce  accurate  representations  of  the  true  system  measurements. 
This  model  will  be  called  the  true  system  model  and  denoted  with  the  subscript  T.  Thus  the 
representations  of  the  true  system  measurements  will  be  denoted  Zr(/,),  and  we  are  assuming  that 
these  true  system  model  representations  portray  the  true  system  measurements  so  that  z(/,)  =  Zj-Ct,). 
Now  assume  that  the  true  system  model  is  a  discrete-time  equivalent  system  model  [33]  of  the  form 

^ r ^i-i) * ^ 
z('/)  =  z^it^)  = 


where 


Xf  is  the  true  system  state  vector 

is  the  slate  transition  matrix,  the  discrete-time  equivalent  of  the  true  system  dynamics 
matrix 

Br  is  the  discrete-time  equivalent  of  the  true  system  control  input  matrix 
«  is  the  system  input  vector 

Gj-  is  the  discrete-time  equivalent  true  noise  input  matrix 

Wj-  is  an  additive  white  discrete-time  dynamics  noise  input  with  zero  mean  and 


0,  ,.*tj 


(2) 
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Zf  is  the  true  system  measurement  vector 
Hj  is  the  true  system  output  matrix 

Vj.  is  an  additive  white  measurement  noise  input,  independent  of  Wj-,  with  zero  mean  and 


t.  *  t. 


(3) 


We  are  assuming  that  the  sampling  rate  is  sufficiently  high  so  that  the  system  modes  and 
noise  models  can  be  considered  constant  over  several  sampling  periods.  Thus  the  Bj,  Gj-,  Hj-, 
Qj,  and  matrices  are  assumed  to  be  time  invariant  over  the  duration  of  the  simulations  used  in  this 
research.  Also,  we  have  assumed  that  the  system  control  input  is  held  constant  over  each  time 
sample. 

The  assumption  that  the  system  model  is  time  invariant  allows  us  to  consider  using  a  steady 
state  Kalman  filter  model  (and  eventually  a  steady  state  constant-gain  Kalman  filter  for 
implementation),  which  will  be  denoted  with  the  subscript  k.  Thus  we  have 


where  x*  is  the  Kalman  filter  model  state  vector 

is  the  Kalman  filter  model  state  transition  matrix 
is  the  Kalman  filter  model  control  input  matrix 
u  is  the  system  input  vector 

is  the  Kalman  filter  model  noise  input  matrix 
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is  an  additive  white  discrete-time  dynamics  noise  input  used  in  the  Kalman  filter  model, 


with  zero  mean  and 


e*.  ^ 

0.  /,  ^  t. 


(5) 


is  the  Kalman  filter  model  measurement  vector 
is  the  Kalman  filter  model  output  matrix 

v*  is  an  additive  white  measurement  noise  input  that  is  used  in  the  Kalman  filter  model.  This 
noise  input  is  assumed  to  be  independent  of  w*,  and  zero-mean  with 


R. 


'i  = 


t,  *  t, 


(6) 


Note  that  the  Kalman  filter  model  and  the  truth  model  are  both  linear  models,  but  the  dimensionality 
of  these  two  models  may  not  necessarily  be  the  same.  In  most  cases,  the  Kalman  filter  model  is  a 
reduced  order  (the  number  of  Kalman  filter  model  states  is  often  a  subset  of  the  unth  model  states) 
version  of  the  truth  model. 

The  Kalman  filter  algorithm  uses  this  model  to  define  time  propagation  and  measurement 
update  equations  of  the  Kalman  filter  state  estimates  and  the  Kalman  filter  state  estimate  covariance 
matrix.  The  Kalman  filter  state  estimate  propagation  equation  based  on  the  Kalman  filter  model  is: 

-^kCV)  (7) 
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where  is  the  Kalman  filter  state  estimate  vector 


Zhitf)  is  the  Kalman  filter  estimate  of  the  measurement  vector  before  it  becomes  available 
ti  is  the  time  just  before  the  measurement  update  at  the  ith  time  sample,  and 
r, .  /  is  the  time  Just  after  the  measurement  update  at  the  (/  -  1)  time  sample, 
and  the  state  estimate  covariance  matrix  propagation  equation: 

PkC^i)  -  Gj  (8) 

The  Kalman  filter  state  estimates  are  updated  using: 

(9) 


where  the  Kalman  filter  gain  is: 

and  the  Kalman  filter-computed  residual  covariance  matrix  A*  is: 


where  is  the  measurement  noise  covariance  matrix  used  in  the  Kalman  filter  model. 

The  Kalman  filter  residual  vector,  shown  in  Eq  (9),  is  defined  as: 

which  is  simply  the  difference  between  the  measurements  (z)  and  the  Kalman  filter  estimates,  based 
on  its  model,  of  those  measurements  before  they  are  taken  x^{t.') ). 

The  Kalman  filter  state  estimate  covariance  matrix  is  updated  using: 
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Pk(tn  =  Pk(*n-^^kOi)B,p,(tr). 


(13) 


Poor  numerical  characteristics  of  these  equations  (particularly  the  last  one)  require  that  this  algorithm 
be  implemented  using  the  U-D  (unitary  upper-triangular  -  diagonal)  covariance  factorization  form. 
The  development  of  this  factorization  form  can  be  found  in  Maybeck  [33]. 

The  steady  state  values  of  the  Kalman  filter  estimate  of  the  state  covariance  matrix  can  be 
precomputed  by  iterating  Eqs  (8),  (10),  (11),  and  (13)  until  steady  state  of  the  covariance  and  gain 
matrices  is  reached.  Once  this  value  for  the  state  covariance  matrix  is  found,  the  steady  state  Kalman 
filter  gain  and  die  steady  state  Kalman  filter  residual  covariance  matrix  are  computed  using  Eq 
(10)  and  Eq  ( 1 1 ).  With  this  steady  state  implementation,  the  state  covariance  matrix,  the  steady  state 
Kalman  filter  gain,  and  the  steady  state  Kalman  filter  residual  covariance  matrices  are  assumed  to  be 
constant  and  therefore  do  not  need  to  be  computed  in  real  time.  The  steady  state  Kalman  filter 
equations  become: 

•^*(t, )  =  ^*•’^*(*1-1  )*  ^*“('m)  (14) 

for  propagating  die  state  estimates  and 

for  updating  the  state  estimates. 

3.2.2  Nomenclature  for  Representing  Mismodeling.  We  will  be 

investigating  the  effects  of  an  incorrect  Kalman  filter  model  on  the  filter's  residual:  Specifically  we 
will  look  at  incorrect  modeling  of  the  state  transition  matrix,  $,  the  output  matrix,  H,  and  the  input 
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matrix,  B.  We  also  are  assuming  that  the  Kalman  filter  model  and  the  truth  model  dynamics  noise 
strength,  Q,  measurement  noise  strength,  R,  and  noise  input  matrices,  G,  are  equivalent.  These 
conditions  were  chosen  because  they  commonly  occur  in  failure  detection  applications  where  MMAE 
is  used. 

One  such  application  is  identifying  flight  control  failures,  where  a  mismodeled  input  matrix 
is  used  to  represent  an  actuator  failure  and  a  mismodeled  output  matrix  would  be  a  sensor  failure. 

For  example,  a  single  actuator  failure  can  be  represented  by  zeroing  a  column  of  the  input  matrix,  B. 
The  result  is  that  the  element  of  the  system  input,  u,  that  corresponds  to  that  column  of  B  will  have 
no  effect  on  the  system  dynamics,  which  is  exaclty  what  we  would  expect  for  a  failed  actuator. 
Likewise,  a  single  sensor  failure  can  be  represented  by  zeroing  a  row  of  the  output  matrix,  H.  The 
result  is  that  the  element  of  the  measurement  that  corresponds  to  the  zeroed  row  of  the  output  matrix 
will  not  consist  of  a  linear  combination  of  the  system  states,  but  will  only  have  the  additive  white 
measurement  noise  v^,  which  is  precisely  the  expected  result  of  a  failed  sensor. 

For  example,  assume  that  the  truth  model  represents  the  case  of  a  single  actuator  failure  and 
the  Kalman  filter  model  assumes  a  non-failed  system;  the  obvious  result  will  be  a  difference  between 
the  truth  model  and  the  Kalman  filter  model.  In  this  case  the  difference  occurs  in  the  column  of  the 
system  input  matrix  that  corresponds  to  the  failed  actuator.  The  truth  model  would  have  a  column  of 
zeros,  while  the  filter  model  would  have  the  non-failed  terms  (some  would  be  nonzero).  Similarly,  a 
failed  sensor  truth  model  would  result  in  a  difference  in  the  row  of  the  output  matrix. 

A  mismodeling  of  the  state  transition  matrix  can  occur  when  the  filter  model  is  based  on  a 
certain  operating  point  (altitude,  velocity,  center  of  gravity,  etc.),  while  the  true  operating  point  is 
different.  The  system  dynamics  are  directly  related  to  the  operating  point.  Thus,  there  would  be  a 
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difference  between  the  filter  model  and  truth  model,  certainly  in  the  state  transition  matrix  and 
possibly  in  the  system  input  matrix. 

In  most  applications  of  Kalman  filtering,  a  designer  will  decrease  the  order  (the  number  of 
states)  of  the  Kalman  filter  design  model  from  that  of  the  truth  model,  to  reduce  the  computation 
loading  required  to  implement  the  Kalman  filter.  This  is  done  by  eliminating  certain  states  that  do 
not  significantly  impact  the  fidelity  of  the  filter  model,  and  then  appropriately  adding  pseudonoise  to 
the  filter  model  to  account  for  the  added  uncertainty  caused  by  the  elimination  of  these  states.  The 
result  of  eliminating  these  states  is  an  incorrectly  modeled  state  transition  matrix,  output  matrix,  and 
input  matrix.  For  example,  if  we  were  eliminating  the  n“*  state  in  the  following  model: 


W) 

‘t’l.l-'t'l.n-l 

= 

*^71-1.1"  *t’n-l.B-l  ‘t’n-l.fi 

^l(^l') 

*l,r"''l.B-l  *1.« 

^«-l'(h) 

= 

1,  r"  l.n- 1 

^B».  1  "*  ^Bi,  B  - 1  n 

■\r-, 

hr 

B-1.  I  ■' 

'■  *B  -1.  r-1 

K-y.r 

'  ^B.  r- 1 

hr 

“i  (^1-1 ) 


) 


(16) 


we  would  need  to  zero  out  the  n'"  column  and  n”"  row  of  the  state  transition  matrix,  the  n*  row  of  the 
system  input  matrix  and  the  n'*'  column  of  the  output  matrix.  This  results  in: 


^l(V) 

't’l.l-'f’l.B-l 

0 

^»-l(h’  ) 

*^B  -  1,  1  *^*B  -  1,  B  -  1 

0 

^B-l  ('<-1 

) 

0 

0-0 

0 

^b('i-i)J 

^'l  ( ) 

*l,r"^l.n-l 

0 

^i(6) 

^Bl  -  1  (  ) 

= 

^rn-l.l  ■■■  ^Bi-1,B-1 

0 

(  h'  ) 

^Bt,  I  "  *B^  B  -  1 

0 

which  essentially  eliminates  the  n"*  state. 
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This  method  of  eliminating  states  results  in  partitions  of  zeros  in  the  model,  which  most 
designers  would  simply  drop  from  the  model  so  that  the  truth  model  and  the  Kalman  filter  models  are 
of  differing  dimensions.  Additionally,  in  many  cases  the  designer  will  aggregate  states  so  that  there 
will  not  be  a  one-to-one  correspondence  between  the  states  of  the  truth  model  and  the  Kalman  filter 
models.  An  elegant  method  that  allows  this  state  reduction  is  described  by  Sheldon  [55,  56].  He 
uses  a  transformation  matrix  that  actually  removes  or  combines  various  truth  model  states,  so  that 
the  truth  model  and  the  Kalman  filter  model  differ  in  dimensionality.  To  simplify  the  development  of 
this  research,  I  have  chosen  to  keep  the  dimensionality  of  the  truth  model  and  the  Kalman  filter  model 
the  same.  Further  research  is  needed  to  develop  the  hypothesis  testing  techniques  using  the 
transformation  matrix  method  described  by  Sheldon. 

We  introduce  the  following  definitions: 

^B^  &  By-  -  B^  =  By-  IHB^ 

AH,  A  By-H,  ^  B,^By-  AH, 

A^,  A  <t>y-  <l>,  -  <S>,^  9y-  ACf,  (Jg) 

where  we  are  implicitly  assuming  that  the  dimensions  of  the  filter  model  and  the  truth  model  are  the 
same. 


Using  this  definition  for  the  examples  above,  we  see  that  in  most  cases  rather  sparse 
matrices  result.  For  instance,  for  the  actuator  failure  where  the  difference  between  the  truth  model 
and  filter  model  is  due  to  a  zeroed  column  of  the  B  matrix,  the  result  would  be; 
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Similarly,  for  a  sensor  failure  where  the  difference  between  the  truth  model  and  the  filter  model  is  due 
to  a  zeroed  row  of  the  H  matrix,  the  result  would  be: 
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Initially,  we  will  assume  that  the  true  system  model  is  accurately  represented  in  one  of  the 
Kalman  filters  of  the  MMAE.  This  assumption  allows  us  to  define  the  true  system  Kalman  filter 
gain,  which  is  the  steady  state  Kalman  filter  gain  of  the  filter  that  uses  the  true  system  model. 

The  Kalman  filter  that  uses  the  true  system  model  will  be  called  the  true  filter.  Later,  we  will  relax 
this  assumption  and  discuss  the  effect  on  the  performance  of  the  MMAE. 

Equations  (5),  (7),  (8),  and  (10)  show  that  the  steady  state  Kalman  filter  gain  is  a  function  of 
$4,  G*,  W*,  Q*,  and  /?*.  Note  that  it  is  not  a  function  of  thus  any  mismodeling  in  the  input  matrix 
will  not  change  the  steady  state  Kalman  filter  gain.  However,  any  difference  between  and  Hj,  or 
<^4  and  will  result  in  a  difference  in  Kalman  filter  gains  between  the  true  filter  and  the  filter  with 
the  mismodeling.  Therefore  we  define 

ATj  k  Kj  -  Kj.-  (21) 

We  now  define  the  error  between  the  Kalman  filter  state  estimate  and  the  true  state  as: 

€4(1:)  A  xj.(r,)  -  :«*(/,■ ),  (22) 
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and 


(23) 


For  the  case  in  which  the  Kalman  filter  model  and  the  true  system  model  match,  that  is  the 
true  filter,  the  state  estimate  error  is  defined  as: 


^7"( 


(24) 


and 


®r('i’)  “  *r(h)  - 


(25) 


Note  that  until  the  filter  state  estimates  are  updated  using  z  (/,),  the  best  estimate  of  jrj.(f,)  is 
Xj-itf).  Therefore: 

•*r(h)  }  =  ■^r(h  ) 


(26) 


where  E  {  •  }  is  the  conditional  expectation  operator,  which  is  conditioned  on  Z(t,.. , ),  the 
history  of  measurements  up  to  and  including  time 

The  true  filter  state  estimate  covariance  matrix  is  defined  as: 

f’r(*i  )  “  ^ z(t,  [•*r(*i)  ■  )][*r(h)  ‘  ■^r(h  )]  } 

"  ^r('/  )  ®r('i  )  }• 


(27) 


3.2.3  Covariance  of  the  Residual.  We  will  derive  the  covariance  of  the  residual  by 
first  deriving  an  expression  for  the  mean  of  the  residual.  Starting  with  Eq  (9)  we  get: 
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(28) 


{  '■*(*,)  } 


"  }  ■  ^Z(.t,.i){  ^k  ^k(^i  )  }  ■ 


Note  that  we  are  conditioning  the  expectation  of  the  measurement  history  up  to  the  r, . ,  time  sample. 
Eqs  (4),  (6),  (7),  (8),  and  (9)  show  that  the  state  estimates  in  the  second  term  of  Eq  (28)  are  simply  a 
function  of  the  Kalman  filter  model  and  the  previous  measurements.  Therefore,  if  we  are  given  this 
measurement  history,  the  state  estimates  are  not  random  and  can  be  directly  computed,  so  the  only 
random  variable  in  Eq  (28)  is  z(t,).  Thus  we  have; 

'■*('■)  }  =  2('i)  }  -  ^k(‘i  )•  (29) 

Now  we  use  Eq  (1)  to  get: 

^ ^ ^r*r(^<)  *  ’'(^)  } 

=  {  Xj.(t^)  I  .  v(/()  I 

=  fl^r^r(V)  *  0 

^r('<  )•  (30) 

Now  to  compute  the  conditional  covariance  matrix  of  the  residual,  conditioned  on  the 
measurement  history  up  to  the  f . time  sample; 
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-  -ff*  ^i(V)  ^Z(^,,){z('.)'’}  ^  B,  £,itr)  £,(t:f  Hj 

■  [^Z(^.,){^^  }  '  ^Z(^.i){^  (  ^■)} 

-  B,i,(tr)  ^z(^.,){z ('.)’■ )  *  B,  £,it-)  £,it-f  ff/j 

=  ■  ■®Z((l,,){^^'«'^}  ^Z(<i.,){z('<)  } 

"  '""'’Zd,.!)  (32) 

This  last  expression  tells  us  that  the  covariance  matrix  of  the  residual  is  n^  dependent  on  the  Kalman 
filter  model,  it  is  only  dependent  on  the  covariance  of  the  measurements.  We  would  expect  this 
because,  as  noted  earlier,  the  only  random  part  of  the  residual,  given  the  measurement  history  up  to 
t,-. ,,  is  z(t, ).  Thus,  in  the  MMAE  filter  bank,  all  of  the  Kalman  filter  residuals  would  have  the  same 
covariance  since  they  all  use  the  same  measurement  vector,  z(/,),  and  measurement  history,  Z(t 

To  get  a  more  explicit  expression  for  the  covariance  matrix  of  the  residual,  we  compute  the 
conditional  covariance  matrix  of  z(t.)  and  use  Eq  (1)  to  get: 

^  zUi.oi  \  Bt^t^  ‘i )  ’  ''r^  *1^]  [Bj-^-p(.ti)  ]  -  Bj.£.j.it^  )  £j.(t^  )  HjF 

=  Bj.E  jffj.  ♦  j.(r^  ) 

*  ^Z(^_,){'’r(  ^  )}  )  B-p  *  ,)  {  ) 

-  Hp£pit:)  £pit:fHp^ 

~  B pE Xp{t^)  ^Hp  *  Rpit.)  -  Hp£p{t.  )  £p{t.  Hp^ 
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-  £^itr)  £^(t-f  .  £j.(t;)  £j.(trf  ]h/  .  Rj.(t.) 

=  ^tI  ^Z(<i.,){ 

-  ;ej.(  V)  "  ^r(V)  ^r(V)^  ]^/  *  ) 

"  ^r[  ^z(^.,){*r(^)  '  ^z(^_,){  *r(^)  }  ) 

- )  ^z(<i_,){  *r( ^  )  ^r( *i  )  j'^r 

+  i?  j-(  ) 

=  BtE  2(t,_^){\^T^^i^  ~  ^\  \  ■  ■^r(  ^  )  ]  )^T  * 

=  nj.Pj.(t;)H/ 1 

Z(<i.,)(''*^'<)}  °  (34) 

where  Aj.(f,)  denotes  the  Kalman  filter  residual  covariance  matrix  at  time  /,.  We  now  assume  that  the 
Kalman  filter  has  reached  steady  state  so  that  the  covariance  matrix  is  constant.  We  denote  the 
covariance  matrix  of  the  residual  of  the  steady  state  Kalman  filter  that  is  using  the  true  system  model 
as  A-p.  This  result  tells  us  that  the  steady  state  covariance  matrix  of  any  Kalman  filter  residual  is 
independent  of  the  Kalman  filter  model  and  can  be  precomputed  simply  by  calculating  the  steady 
state  residual  covariance  of  a  Kalman  filter  with  the  true  model. 

3.2.4  Mean  of  the  Residual.  Earlier  it  was  shown  that  the  conditional  covariance  of 

the  residual  from  any  of  the  Kalman  filters  in  the  MMAE  is  dependent  only  on  the  conditional 
covariance  of  the  actual  measurements;  thus  any  mismodeling  in  the  design  model  upon  which  the 
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filter  is  based,  has  no  effect  on  the  conditional  covariance  of  the  residual.  We  will  now  examine  the 
effects  of  mismodeling  on  the  conditional  mean  of  the  residual. 

First,  we  expand  one  of  the  terms  in  Eq  (12)  using  die  definition  in  Eq  (1)  to  get: 

'■*(^)  ^  *  M^)]  - 

.  [ffj.Br-  [Bt-  ^B^)[Bj.-  AB*)]ii(r,._,)  *  ffj.Gj.WjCt._^)  .  v^Cr,) 

=  ffj.9j.[xj.(t^_^)  -  •  [ffr^^k  *  ^B^^j  -  A  A 

♦  (.ffrABj  *  -  Afl’j  A  .  ffj.Gj.w^Cti.i)  * 

=  ^  ^ **)^*( '/-i ) 

,  [ffy^Bk  *  ^B^Bj.  -  Afl’*  A<I>^)h  (r^., )  .  ffj.Gj.w^(t..^)  .  Vj.(q) 

(35) 

We  take  the  expectation  of  this  to  find  the  mean  of  the  residual; 

♦  (fl’j-ABj  •  -  Afl-jA^j)!/  (/,.,  )  .  Hj.Gj.Wj(t._^)  .  Vr(/,)} 

.  ’  ^B^Bj. 

.  AF^B^-  AB^A<t,)ii(r..,)  (36) 
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To  derive  an  expression  for  the  first  term,  we  take  the  conditional  expectation  of  Eq  (23)  and  then 
step  back  one  data  sample  in  time: 


"  ^zct,. 

i)f 

*7-(*,) 

,.l 

^7"(  )  ■ 

[£,Ct;)-  K, 

:Z(*l)]} 

=  ■®'z(^. 

,.l 

J^T-C  *i)  - 

{l-K,H,)£,(t;)-K,[H^x 

r(*/)  *  '’r(*i)]} 

"  -^zcii. 

.)! 

T)xT(t,)-[I-K,H,)£,Ct;) 

}  -  ^k^Z(.t,, 

.,){^r(*i)} 

°  ^Z(^. 

r)(  *1-1 )  *  -^r"  ( *1-1 ) 

*  ^r’*’rf(*i- 

i)] 

- 

*t^t(*;-i)^^***(*i-i)l} 

■C* 

N 

II 

i)( 

[i-K,n 

r)«>rM*.-i)- 

)} 

*  (f  -  *1-1 ) } 


Now  we  substitute  in  the  definitions  from  liq  (18): 


=  (^■^*^r)*r^z((i,,){*r(*i-i)  ‘  ^t(*i-i )} 

.  [{I-K,H^)B^-  {l-K,Hj.)Bj..  {l-K,H^)liB,-  K,^^,B^]u(t,_^ 
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(38) 


We  will  assume  that  until  a  certain  point  in  time  the  Kalman  filter  model  matched  the  true 
system  model.  For  failure  identification  applications,  this  point  in  time  would  be  the  failure  time,  tj. 
Up  until  that  time,  all  of  the  modeling  errors  are  zero  and 


(  h'  ) 


V/, 


<y 


(39) 


Therefore  the  mean  of  the  state  estimate  error  will  also  be  zero.  This  assumption  gives  us: 


0 

,  for  r,<=r^ 

.  for 

(40) 

Thus,  to  compute  the  mean  of  the  residual,  we  keep  a  running  calculation  of  the  mean  state 
estimate  error  using  Eq  (40) ,  and  then  compute  the  residual  mean  using: 
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We  can  now  use  this  equation  to  look  at  the  mean  of  the  Kalman  filter  residual  for  several 


special  cases  of  mismodeling. 

3. 2. 4.1  No  Mismodeling,  We  will  first  look  at  the  case  where  the  Kalman 


filter  model  and  the  true  system  model  match.  For  this  case  we  get: 


A£*=0 

0 

A$*=  0. 


(42) 


Using  Eq  (40)  and  Eq  (41 )  we  find  that: 

"  =  0  t  (43) 

which  Maybeck  [33]  has  also  shown  using  a  different  development. 

3. 2. 4. 2  Mismodeled  Input  Matrix.  A  mismodeled  input  matrix,  which  is 
used  to  model  an  actuator  failure  in  flight  control  failure  identification  applications,  would  result  in: 


Bt*  B, 

Ht  = 


ABj  *  0 

AF*=  0 
A3.,=  0. 


(44) 


Using  Eq  (40)  we  get: 

)}  = 


0  ,  for  t.  s 

[I-K^Hr)^B^u(tf)  ,  for  = 


(45) 
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and  from  Eq  (41)  we  get: 


(46) 


3. 2.4. 3  Mismoddled  Output  Matrix,  a  mismodeled  output  matrix,  which 
is  used  to  model  a  sensor  failure  in  flight  control  applications,  would  result  in: 


Bt 

Hr 


B 

H 

i 


0 

0 

0. 


(47) 


Using  Eq  (40)  again,  we  get: 


0 

-  -  K^AH^B^u  (t^) 

-  AH  1^  ij  _  ^)  -  Ki^AHi^Bi^u(_t.^) 


,  for  i 

>  f°’-  ',  =  '/•! 


for  , 

(48) 


which  can  be  rewritten  as: 


0  ,  for  t.  s  ty 

,  for 

-K,AH,£,(t,) 


(49) 


Using  Eq  (41)  we  get: 
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=  ,){  ^J[(  *1-1  )} 


®  (  *1-  1  ) 


Afl'*5r“(*.-i) 

(50) 


3. 2.4.4  Mismodeled  State  Transition  Matrix,  a  mismodeied  state 

transition  matrix,  which  can  be  used  to  model  an  imperfect  modeling  of  the  aircraft  modes  in  flight 
control  fuiluic  identification  applications,  would  result  in: 


Bt 

Br 


AB,  =  0 

A^*=  0 

A4>j  #  0. 


(51) 


Using  Eq  (40)  again,  we  get: 


^Z(f,.,){®*(*i  )}  = 


0  ,  for  t.  s 

,  for 


(52) 


Using  Eq  (41)  we  get: 

"  ^r^r^z((i_,){^*(*i-i )}  *  *t^*(  *<-i  )•  (53) 

3.2.5  Kalman  Filter  Bank  Implemented  Using  a  Single  Residual.  Indus 

section  we  propose  a  linear  transform  that  operates  on  the  residual  from  one  Kalman  filter  and 
produces  the  equivalent  residual  of  another  Kalman  filter.  The  proposed  structure  is  shown  in  Figure 
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Figure  5  -  Alternative  Method  of  Computing  a  Second  Kalman  Filter  Residual 


5.  A  single  Kalman  filter,  denoted  fdter  k,  is  implemented  and  its  residual  and  state  estimates,  along 
with  the  control  input  and  measurement  vectors,  are  used  by  a  linear  transform  to  produce  the 
equivalent  of  the  residual  from  another  Kalman  filter,  denoted  filter  j.  This  transform  only  requires 
knowing  the  difference  between  the  filter  models  and  docs  not  require  the  implementation  of  filter  j. 
This  structure  can  be  extended  to  form  the  equivalent  of  the  Kalman  filter  bank  shown  in  Figure  4. 
This  alternative  structure  is  shown  in  Figure  6,  where  a  single  residual  is  used  to  form  the  equivalent 
residual  from  all  the  other  Kalman  filters. 

This  structure  could  produce  significant  computational  savings  for  many  MMAE 
applications.  We  illustrate  this,  we  will  compare  the  operations  count  (number  of  multiplies  and 
additions)  that  are  need  to  implement  the  usual  Kalman  filter  bank  structure  with  the  count  needed  to 
implement  this  structure.  We  will  not  assume  a  specific  structure  (i.e.,  diagonal,  block  diagonal, 
partitioned,  etc.)  of  the  matrices  used  to  implement  these  structures,  except  for  the  AB,  AH,  and  A<^ 
matrices.  Previously,  we  showed  that  these  matrices  are  sparse,  usually  with  only  one  nonzero  row  or 
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Figure  6.  Alternative  Multiple  Model  Adaptive  Estimation  Filter  Bank  Structure 

Using  Equivalent  Residuals 


column.  We  will  develop  the  operations  count  for  the  fully  implemented  Kalman  filter  bank  in  this 
section,  and  compare  this  to  the  operations  count  for  the  case  of  a  mismodeled  system  input  matrix 
(simulating  an  actuator  failure)  in  Section  3.2.5. 3  and  a  mismodeled  output  matrix  (simulating  a 
sensor  failure)  in  Section  3.2.5.4. 

We  start  by  defining  the  various  matrix  dimensions  as  follows:  -  number  of  states, - 
number  of  system  inputs,  and  -  number  of  measurements.  The  steady-state  Kalman  filter  state 
estimates  are  updated  using  Eq  (15),  where  the  dimension  of  is  by  n„,  the  state  estimate  vector 
is  n^,  and  residual  vector  is  n^.  Thus,  we  will  have  multiplies,  and  (n„-l)+«^  = 
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additions.  The  state  estimates  are  propagated  using  Iiq  (14),  where  the  dimension  of  O*  is  by  n^, 
is  by  n,,  and  the  system  input  vector  is  Eq  (14)  requires  «,  +  multiplies,  and  {n^  - 
l)+«^  (n^  -l)+n^  =  additions.  Lastly,  Eq  (12)  is  used  to  form  the  current  time  sample  of  the 
residual,  which  requires  multiplies,  and  {n^  -!)+«„  -n^  additions.  Thus,  the  total 
operations  coimt  is  (n^  +  «,  +  2n„  )  multiplies  and  (n^  +  n,  +  2n^  -1).  The  dimension  for  the 
specific  application  that  was  implemented  for  this  research  were  =19,  n,  =6,  and  =  9,  which 
gives  us  an  operations  count  of  817  multiplies  and  798  additions,  per  elemental  filter.  This 
operations  count  will  be  compared  to  the  count  using  this  alternate  structure  in  the  following 
sections.  Now  we  will  compute  the  linear  transform. 

3.2.5. 1  Model  Difference  Nomenclature.  This  development  will  be 

similar  to  Section  3.2.2,  except  that  we  will  be  developing  the  difference  between  two  Kalman  filter 
models  instead  of  between  the  true  system  model  and  a  Kalman  filter  model.  We  will  use  the 
subscripts  j  and  k  to  denote  the  two  different  Kalman  filter  models. 

We  are  examining  model  differences  in  the  state  transition  matrix,  O,  the  output  matrix,  H,  and 
the  input  matrix,  B.  Therefore,  wc  introduce  the  following  definitions; 


A 

Bu 

-B.  ^ 

5.- 

Bu- 

A 

-H.  ^ 

H.- 

B,- 

A 

-  - 

$.= 

A 

-K.  ^ 

Jf.= 

In  Section  3.2.2  we  defined  the  error  between  the  Kalman  filter  state  estimates  and  the  true 
state  using  the  nomenclature  e*.  In  Section  3.2.4  we  showed  that  the  mean  error  between  the  Kalman 
filter  state  estimates  is  important  in  computing  the  mean  of  the  residual  and  must  be  computed  at 
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each  time  sample  to  compute  the  mean  of  the  residual.  In  a  similar  manner,  we  define  the  difference 
in  state  estimate  errors,  which  is  algebraically  equivalent  to  the  difference  in  state  estimates  as: 


=  [  ] 

-  =  Ac.,(/;).  (55) 

3. 2.5.2  Computation  of  the  Residuals.  Our  goal  is  to  replace  the >-th 

Kalman  filter,  in  the  MMAE  filter  bank,  using  the  residual  from  the  ^-th  Kalman  filter.  We  assume 
that  the  Kalman  filter  states  from  filter  k  and  the  residual  from  filter  k  are  computed  at  each  time 
sample  and  are  available  for  computations.  Also  the  control  input  and  the  measurement  vectors  are 
available  for  computations.  We  assume  that  the  model  differences  are  only  between  the  input 
matrices,  output  matrices,  and  state  transition  matrices  (or  some  subset  thereoO-  Thus,  we  desire  to 
derive  an  expression  for  the  residual  of  filter  j,  in  terms  of  the  state  estimates  of  filter  k,  the 
difference  between  the  Kalman  filter  state  estimates  of  filter  k  and  j,  and  the  known  model  matrices. 

We  start  with  the  definition  of  the  residual  for  filter /: 

rjit,)  & 

=  z(/,)-  (/,.,)] 

=  z(/,)  -  -  £r.B. «(/,.,) 

=  z(/,)  - 

-  Z(^)  -  .  ®.  A  e.*( .  (5^) 

This  expression  can  be  used  to  compute  the  equivalent  residual  from  filter  j  using  the  state 
estimates  from  filter  k,  the  difference  between  the  state  estimates  of  filter  k  and  filter  j,  and  the 
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control  input  and  measurement  vectors.  We  will  now  derive  an  equivalent  expression  that  uses  the 


residual  from  fdter  k. 


r.(r,)  =  z(r,)  - 

=  z(/,)  - 

=  Z(q)  - 

Thus,  we  get: 

'•y(<,)  =  '•*(*,)  *  Bj4fjAej^it._{) 

This  may  seem  to  be  a  more  complex  expression  than  Eq  (56),  but  later  we  will  see  that  in  most 
cases  this  simplifies  dramatically. 

To  compute  the  difference  in  state  estimates  between  filter  k  and  filter  j,  we  can  compute 
each  of  these  estimates  and  then  difference  them.  However,  by  doing  so  we  would  need  to  implement 
filter  j,  which  would  defeat  the  purpose  of  this  derivation  since  we  would  have  the  residual  directly 
from  that  filter.  When  the  MMAE  is  initialized,  the  difference  between  the  Kalman  filter  state 
estimates  is  either  zero  or  known.  Therefore,  we  need  to  derive  a  recursive  expression  for  the 
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difference  in  state  estimates  between  filter  k  and  filter y  ( A  ( r/ ) ).  We  could  then  compute 
A  )  at  each  time  sample  and  then  the  desired  residual  for  the  filter  j. 


We  start  with  Eq  (55)  and  simply  decrement  back  one  time  sample: 


=  {  \  )  -  {  A. it:)  .  K.  [  z(l,)  -  H.A.it;)]  } 

=  (  /  -  K,n^  )A,it:)  -{I-  K.H.  )Aj(t[)  .{K,-  Kj  )z(/,.) 

=  (7  -  K,H,  )^,A,(t,_,')  -  (7  -  KjHj  )^jA. 

M  (^  -  ^jBj  ]«(V,)  •  (^*  -  )z('.) 


Now  we  substitute  in  the  definitions  from  Eq  (54)  and  Eq  (55): 

Ae^.,(l/)  =  (7  -  K,H,)^,A,(t,,{)  -  {I  .  KjHj)<S>.\A,(t,_{)-iiej, (,,_;)] 

M  (  ^  -  ^kB,  )B,-  (I-  KjHj  )Bj  ]  «  (  1,. .  )  ,  Ar,  .Z  (  ,,  ) 

=  (7-ir,ff,)«,Ae,,(t,.,-) 

M  (^  -  K  -  (^  -  )^J  ]“('.-!>  ^  (60) 

which  is  the  desired  recursive  relationship. 

Eq  (58)  and  Eq  (60)  certainly  appear  to  be  very  complex  expressions,  but  under  most  model 
differences  these  expressions  will  simplify  dramatically.  We  now  look  at  some  specific  model 
differences  to  demonstrate  this  simplification. 
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3. 2. 5. 3  Different  Input  Matrix  Models.  We  now  assume  that  the  two 

Kalman  filter  models  differ  only  in  the  input  matrix.  This  would  represent  different  actuator  failure 
conditions  for  the  filter  control  failure  identification  application.  This  gives  us; 


B,  *  Bj 


0 

=  0 
AJT,.  =  0 

AB,j  ^  0 


(61) 


Under  these  assumptions  Eq  (60)  becomes: 

Ae.*(/,-)  =  {T  - 

.[{l.KjHj)B,.{l-  K.H.  )Bj  ]«(/,.,)>  AK.jZ  ( t, ) 

Ae.*(t/)  -{I  -  KjHj  )<1>,A€.*(1,.,-) 

*  -  KjHj)AB^jU(t,_j)  (52) 


and  Eq  (58)  becomes: 

'■yC'i)  =  *  B.AB^jU{t._^).  (53) 

Returning  to  our  flight  control  failure  detection  example,  we  find  that  we  could  construct  the 
equivalent  residuals  of  a  bank  of  Kalman  filters,  each  with  a  different  actuator  failure  model,  by 
using  a  single  Kalman  filter  residual  and  Eqs  (62)  and  (63).  We  will  now  compute  the  operations 
count  and  compare  it  to  the  required  number  of  operations  for  the  fully  implemented  Kalman  filter 
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bank.  We  define  as  the  number  of  columns  of  AB  that  are  not  zero,  thus  a  mismodeling  in  a 
single  column  of  B  will  produce  a  AB  with  only  one  nonzero  column,  and  rig-l.  The 
which  is  by  n^,  and  (i-KjHj)  AB^j  terms  in  Eq  (62)  are  precomputable,  and  therefore  will  not  be 
computed  at  each  time  sample.  The  (I-KjHj)  ab^j  term  is  a  sparse  by  «,  matrix  with  only  Hg 
nonzero  columns.  To  evaluate  the  first  term  of  Eiq  (62)  will  require  multiplications  and  n/«^-l) 
additions.  The  second  term  requires  rig  multiplications  and  n^{ng-\)  additions.  A  similar  analysis 
of  Eq  (63)  shows  that,  the  second  term  requires  multiplications  and  njn^-l)  additions,  the  third 
term  requires  n„  iig  multiplications  and  n„(«fi-l)  additions,  plus  2n^  additions  to  sum  the  first, 
second,  and  third  terms  together.  Thus  the  total  operations  count  is  {n^  +  +  rig)  +  iig 

multiplications,  and  (rig  +  n^)  +  (rig-l)  additions.  This  results  in  560  multiplies  and  541 

additions  for  the  specific  application  that  we  implemented  for  this  research.  Subtracting  this  count 
from  the  total  operations  count,  developed  earlier,  gives  us  the  total  operations  savings  of 
implementing  this  structure  over  the  fully  implemented  Kalman  filter  structure.  This  savings  is 
(n,  +  rij  -  rig  (ri^  +  nj  multiplies  and  the  same  for  the  number  additions,  which  is  257  for  this 
particular  application.  This  savings  would  be  multiplied  by  the  number  of  filters  that  are  modeling 
actuator  failures  in  the  Kalman  filter  bank,  to  produce  the  total  savings  of  implementing  this  structure 
over  fully  implementing  the  Kalman  filter  bank.  For  this  specific  implementation,  the  savings  is 
about  30%  of  the  required  operations. 

3.2.5.4  Different  Output  Matrix  Models.  We  now  assume  that  the  two 

Kalman  filter  models  differ  only  in  the  output  matrix.  This  will  cause  the  two  filters  to  also  have 
different  Kalman  filter  gains,  as  explained  in  Section  3.2.2.  This  gives  us: 
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Under  these  assumptions  Eq  (60)  becomes: 

{I  - 

M  (^  -  )^t  -  (^  -  )^*  ]“(',-i)  * 

*  [{^  -  -  {I  -  I^jBj)]B,u(t^_,)  ,  iiK,jz(t,) 

=  (/-r.F^.)«.,Ae.*(l,.,’),  AJC*.z(/,) 

.[KjHj  -  K,n,  ] 4* .  [  AT.^.  -  K,H,  ]fi*  «(/,.,) 

-  I  r  I'  rr  \  \  c  (  >  •^  »«'•  -<'<'1 

V  -j-j  ; 

^[KjB.  -  K,H,  /,.,*)  .£*«(/,.,)] 

A€.t(V)  =  (/  -  K.Hj  )  4>,  A  e. *(/,_,•)  *  AX^-zCr.) 

and  Eq  (58)  becomes: 

'M)  =  '■*(',)  *  B.i>.L^.,it._{)  ,  ,  ^H^.B,uit._,) 

=  »•*('*)  »  Bj^jAej^(t,_;)  .  ,  £*«(/,.;)] 

rj(ti)  =  r4(r.)  . 
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Returning  to  our  flight  control  failure  detection  example,  we  find  that  we  could  construct  the 
equivalent  residuals  of  a  bank  of  Kalman  filters,  each  with  a  different  sensor  failure  model,  by  using 
a  single  Kalman  filter  residual  and  Eqs  (65)  and  (66).  We  will  first  develop  the  operations  count 
assuming  that  AK^j  is  not  necessarily  sparse,  and  then  develop  the  operations  count  assuming  that 
is  sparse  with  a  small  nonzero  partition.  We  define  as  the  number  of  nonzero  rows  of 
Since  none  of  the  terms  in  Eq  (65)  are  sparse,  the  operations  count  for  this  equation  is  2n/  + 
multiplications  and  2n/+  («„-l)  additions.  The in  Eq  (66)  is  sparse,  so  only 

multiplies  and  additions  are  need  to  evaluate  the  last  term.  Thus,  the  operations  count  for 

Eq  (66)  is  (riff  +  nj  multiplies  and  +  nj{  n^-l)  +  additions.  The  total  count  for 
implementing  the  equivalent  residual  structure  is  (  2n^  +  2n^  +  n„)  multiplies  and  (  2n^  +  2n^  + 
additions.  This  structure  requires  1083  multiplies  and  1072  additions  for  the 
specific  application  that  we  implemented,  only  because  we  assumed  that  AKj^j  is  not  necessarily 
sparse. 

We  observed  that  AK^j  is  actually  sparse  for  the  particular  application  that  was  implemented 
for  this  research.  This  matrix  would  have  a  small  block  partition  that  would  be  nonzero  for  the 
sensor  failures  that  were  modeled.  We  define  the  dimension  of  the  block  as  n^,  and  expanded  the  last 
term  in  Eq  (65)  to  get  When  the  sparse  structure  of  AK^jmA  AH^jWce 

taken  into  account,  this  term  is  a  sparse  matrix  with  nonzero  rows  and  an  by  nonzero  block 
partition.  This  results  in  a  operations  count  of  multiplications  and  (n^  -1)  + 

{n^  -1)(«^  -1)  additions.  The  total  operations  count  for  Eq  (65)  and  Eq  (66),  when  the  sparse 
structures  of  AK^^j  and  AH^^  are  taken  into  account,  is  +  +  n^J  +  1)  +  (2njc  -  1)  multiplies 

and  («j  +  «„  +  +  1)(  -  1)  +  {2nf,  -  1)  («jf  -  1)  +  2n^  additions,  which  gives  us  615  multiplies  and 

594  additions  for  this  particular  implementation  where  tif,  =  1  and  n^  =  5.  Subtracting  this  from  the 
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operations  count  for  the  fully  implemented  Kalman  filter  version  gives  us  a  savings  of  (n„  +  n,  -  n„ 
-  1)  +  rijj-  2n^  multiplies  and  - 1)  +  «;,  +  1  -  Q-tif,  -  l)(njf  - 1)  -  which  is  a  savings 

of  202  multiplies  and  204  additions.  Thus,  the  equivalent  residual  version  of  the  Kalman  filter  bank 
requires  about  25%  fewer  operations  than  the  fully  implemented  Kalman  filter  bank. 

3. 2,5. 5  Different  State  Transition  Matrix  Models.  Finally,  we 

assume  that  the  two  Kalman  filter  models  differ  only  in  the  state  transition  matrix.  This  will  cause 
uie  two  iiiiers  to  iiave  diileiciu  Kaiiauii  fiiici  gains,  as  cxpiamcu  m  Sccuoii  3.2.2.  Tins  gives  us; 


K,  *  K. 


^  0 


AB*.  =  0 


Aff*.  =  0 


6.K^.  *  0 


(67) 


Under  these  assumptions  Eq  (60)  becomes: 

A6.*(i;)  =  (/-X.^.)4..Ae.,(l,..,-) 

)^*]-(Vi)- 

=  (/-jr,ff,)e,Ae, *(/,.,•) 

.[K.Hk-KkB^  ]Bj  «(/,.,).  AJr,.z(t,) 

=  (/-  Ar.ff.)$.Ae.,(l.._,-)-  Ar,.z(t.) 
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Ae.*(V)  =  (/  -  KjHj  .  AK.jzH,) 

=  {I-KjHj)i>.Aej,Ct,_,^).  AK.jZit.) 

.  {I  -  KjH,)A^,j£,(t,_{)  -  AK,jH,lib,£,it,,{)  .  B,uit,_^)] 

=  (/  -  KjH.  )^jAej,(t,,{)  .{I-  K.Hj  )  A  ) 

*  AJr,.[z(/,)-£r,^,(/,-)] 

A6.,(t;)  =  (/  -  K.H^  )[*;Ae.,(/,,,*)  .  A$,,i,(/,.,-)]  .  AZ*.r,(/,.)  (68) 


and  Eq  (58)  becomes: 

'•/',)  =  »•»(*.)  ^  ^j^j^^jkiti-{)  *  Bj^^kj^k(h-\)-  (69) 

We  can  still  construct  the  equivalent  residuals  of  a  bank  of  Kalman  filters,  each  with  a 
different  state  transition  matrix  model,  by  using  a  single  Kalman  filter  residual  and  Eq  (68)  and  Eq 
(69).  The  differences  between  models  is  sometimes  in  only  one  row  and  column  of  the  state 
transition  matrices,  so  A  may  be  a  matrix  with  only  one  nonzero  row  and  column.  When  this  is 
multiplied  with  another  matrix,  the  product  is  no  longer  a  sparse  matrix.  Thus,  these  equations  will 
not,  in  general  produce  a  savings  in  computational  loading.  However,  just  like  the  case  of  a 
mismodeled  output  matrix,  the  equivalent  residual  implementation  may  produce  computational 
savings  for  some  specific  applications.  Thus,  the  computational  savings  must  be  investigated  for 
each  specific  application. 

3.2.6  Residual  Correlation  Kalman  Filter  Bank,  in  this  section  we  propose 

constructing  a  Kalman  filter  bank  with  outputs  that  are  estimates  of  the  power  spectral  density  of 
each  of  the  residuals  from  the  Kalman  filters  in  the  bank.  We  showed  earlier  that,  if  the  Kalman  filter 
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Figure  7.  Residual  from  the  Kalman  filler  with  the  Fully  Functional  Model 
in  the  Presence  of  an  Elevator  Failure. 


model  is  correct,  the  residual  is  a  white  sequence  with  zero  mean,  but  if  the  model  is  incorrect  the 
mean  changes.  This  causes  a  change  in  the  residual  correlation,  but  not  the  covariance  of  the 
residual.  In  Section  3. 2.4. 2,  we  showed  that,  if  the  model  differences  are  in  the  input  matrix,  the 
change  in  mean  of  the  residual  is  a  summation  of  input  terms.  If  the  input  is  a  sinusoid,  then  these 
terms  produce  a  sinusoidal  residual  at  the  same  frequency  as  the  input.  This  effect  is  clearly  evident 
in  Figure  7  where  the  residual  from  the  Kalman  filter  based  on  a  model  that  assumes  a  fully 
functional  aircraft,  shows  the  presence  of  the  elevator  dither  input  when  an  elevator  failure  occurs. 
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Figure  8.  Fourier  Transform  of  the  Pitch  Rate  Residual  from  the  Fully  Functional 
Kalman  FUter  with  Mismodeling  (solid  Une)  and  No  Mismodeling  (dotted  line). 

Since  we  know  the  frequency  of  the  input,  we  can  use  the  spectral  content  of  the  residual  at 
this  particular  frequency  to  indicate  the  presence  of  mismodeling.  Figure  8  is  the  Fourier  transform 
of  200  data  points  of  the  residual  for  the  pitch  rate  residual  element  shown  in  Figure  7  (solid  line), 
along  with  the  same  residual  element  from  the  same  Kalman  filter  when  there  is  no  failure  (dotted 
line).  The  "spike"  in  the  solid  line  occurs  at  the  elevator  dither  input  frequency.  Note  that,  at  this 
particular  frequency,  the  spectral  content  of  the  residual  with  the  mismodeling  is  significantly  greater 
than  the  spectral  content  for  the  correctly  modeled  residual.  This  figure  shows  that  the  spectral 
content  of  the  residual  clearly  indicates  the  presence  or  absence  of  the  mismodeling. 
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3.2.6.1  Basic  Equations.  Kay  [26]  develops  several  spectral  estimation 


techniques;  we  have  chosen  to  use  the  periodogram  [26:  65]  for  this  research  since  its  characteristics 
are  well  researched  and  it  is  used  as  a  conceptual  basis  for  many  other  spectral  estimation  techniques. 
One  version  of  the  periodogram  utilizes  the  fast  Fourier  U'ansform,  which  could  greatly  aid  in 
implementing  this  technique  because  this  transform  has  been  implemented  on  commercial  chip  sets. 
Other  spectral  estimation  techniques  still  need  to  be  researched  to  determine  which  one  or  ones 
produce  the  desired  performance.  The  periodogram  is  based  on  estimating  the  autocorrelation  of  the 
residual  and  then  taking  the  discrete  Fourier  transform  of  the  autocorrelation  to  produce  an  estimate 
of  the  power  spectral  density. 

We  must  first  make  some  assumptions  to  be  able  to  estimate  the  autocorrelation  of  the 
residual.  We  assume  that  the  residual  sequence  is  a  series  of  samples  of  a  stationary  process,  so  that 
the  probability  distributions  do  not  change  with  time.  We  also  assume  that  the  residual  is  ergodic  in 
the  autocorrelation  function,  which  implies  that  the  expected  value  of  the  time-averaged 
autocorrelation  function  is  the  same  irrespective  of  the  length  of  time  averaging. 

For  a  multidimensional  sequence,  the  estimate  of  the  autocorrelation  is: 


AT-l-l/l 


=  ^  E  ^(r,. 

IS 


n  •  0 


The  periodogram  is  the  discrete  Fourier  transform  of  this  sequence,  thus: 


(70) 
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N-\ 

E  A,(l)  cxpi 

N-\  .  W-l-|/| 

=  E  t;  E  'kifi-n-m)  J^^/n 

^  B-o  (71) 


where  N  is  the  number  of  data  samples  that  are  collected  over  time. 

Note  that  by  dividing  by  N  instead  of  (A(  - 1)  or  (A^  -|  /|  -1)  in  Eq  (70),  we  are  using  a  biased 
estimate  of  the  autocorrelation.  We  need  the  MN  factor  to  make  this  an  estimate  of  the  power 
spectral  density. 

If  the  residual  sequence  is  a  scalar  sequence,  than  Eq  (7 1)  can  be  shown  to  be  equivalent  to; 


N 


N-\ 


E 


2 


CXp(  jlTifl) 


(72) 


which  is  simply  the  squared  absolute  value  of  the  N-point  Fourier  transform  of  the  residual  sequence. 
This  is  much  easier  to  implement  and  executes  much  faster  than  Eq  (71),  since  it  exploits  the  fast 
Fourier  transform  routines.  If  the  residual  sequence  is  multidimensional,  we  can  approximate  Eq 
(7 1)  with  Eq  (72)  if  the  cross-correlation  terms  between  elements  of  the  residual  vector  are 
negligible.  Eor  the  specific  flight  control  application  that  we  are  studying,  most  of  these  cross  terms 
are  negligible.  Since  we  are  attempting  to  find  practical  implementations  of  the  MMAE  algorithm, 
and  Eq  (71)  is  computationally  intensive,  we  chose  to  implement  Eq  (72)  for  this  research,  but  we 
will  continue  the  theoretical  development  based  on  Eq  (71),  since  it  is  correct  for  multidimensional 
residual  sequences. 
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Figure  9.  Multiple  Model  Adaptive  Estimation  Algorithm  using  a 
Residual  Correlation  Kalman  Filter  Bank 


We  propose  altering  the  Kalman  filter  bank  by  estimating  the  spectral  content  of  each  of  the 
residuals  using  Eq  (72).  This  structure  is  shown  in  Figure  9  (the  "Hypothesis  Testing  Algorithm" 
block  will  be  specified  in  detail  in  Section  3.3.1).  Kay  [26:  65-66  ]  shows  that  this  can  be  interpreted 
as  filtering  the  residual  with  a  bandpass  filler  centered  at / and  a  3dB  bandwidth  of  l/N,  sampling  the 
output,  and  computing  the  squared  magnitude.  The  IW  factor  is  need  to  make  the  estimate  a  power 
spectral  density.  Clearly,  as  more  data  samples  are  used  (A^  increases)  the  filter  bandwidth  narrows 
and  more  of  the  out-of-bandwidth  noise  is  rejected.  To  use  this  estimate  of  the  power  spectral 
density  in  the  hypothesis  testing  algorithms  that  are  developed  in  Section  3.3,  we  need  the  mean  and 
covariance  of  this  estimate. 
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Note  that  in  Eqs  (71)  and  (72)  we  are  multiplying  two  Gaussian  random  variables  together, 
which  will  yeild  a  chi-squared  distributed  random  variable,  and  then  summing  several  of  the  chi- 
squared  random  variables  together  (we  used  7/ =  100  for  this  research).  Since  each  sample  of  the 
residual  is  independent  and  identically  distributed,  we  can  use  the  Lindeburg-Levy  Theorem  [53]  to 
show  that  the  distribution  of  ❖  converges  in  distribution  to  a  normal  distribution.  We  generated 
histograms  of  using  120  data  points  and  observed  that,  for  the  case  where  the  residual  has  a 
nonzero  mean  (when  a  failure  occurs),  (his  approximation  works  well.  However,  when  the  residual  is 
zero-mean,  then  #  appears  to  be  more  chi-squared  distributed.  We  also  observed  that,  for  the  dither 
input  levels  that  we  used  for  this  research,  the  change  in  mean  effects  of  ^  were  much  greater  than 
the  effect  of  using  an  incorrect  distribution.  In  Section  3.3. 1 . 1 ,  wc  will  approximate  ^  as  Gaussian. 
Further  research  needs  to  be  accomplished  to  characterize  the  distribution  of  t  properly,  for  a  small 
number  of  data  samples. 

3. 2.6. 2  Mean  of  the  Power  Spectral  Density  E.stimate.  To  develop  the  mean 

of  this  estimate,  we  first  need  to  compute  the  autocorrelation  function  of  die  residual.  Earlier,  we 
found  that,  if  there  is  no  mismodcling,  (he  autocnrrcla(ion  ma(rix  is: 


A,(l)  =  i 


^i(O).  for  I  ^  0 


0,  for  1*0 
=  ,4^(0)  6(1) 


(73) 


where  A*(0)  =  and  6(/)  is  the  Kronecker  delta  function  defined  as: 


6(1) 


1,  for  I  ^  0 
0,  for  I  *  0 


(74) 
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We  will  now  allow  the  mismodeling  that  we  studied  in  the  earlier  sections  and  develop  an  expression 
for  the  residual  autocorrelation. 


We  start  by  expanding  the  definition  of  the  residual  (letting  v(r,.)  =  since  the  subscript 
holds  no  information  for  this  application): 


»•*(/,)  = 

=  t.)  4  v(<j. )j  -  ) 

=  Hj.Xj.(t.)  -  4  v(/.) 

=  *  v(/.) 

filial)  -  B  j  )  4  A  (^<)  *  ''(*/) 


We  assume  that  the  mismodeling  has  occurred  over  the  entire  length  of  data  sampling,  all  the  way 
back  to  where  I  is  the  number  of  sample  periods  over  which  wc  are  interested  in  estimating  the 
autocorrelation.  Since  the  residual  is  a  real-valued  sequence,  wc  need  only  estimate  the 
autocorrelation  for  positive  /  because  the  autocorrelation  estimate  for  positive  /  is  the  same  as  for 
negative  /  (A^C/)  =  A^{-1)).  Now  we  use  the  definition  of  the  autocorrelation  of  the  residual,  for  /  ^  0: 

"  )*  )*’'(^)][^*(^-z)  Bj-  *v  ^(r,.;)]  } 

=  h  )^ir(^-z)}^r  *Bj.E2(^f^  )^k  ^B^ 

*  ^z(zi,)  {"*(  h  )  V  ""C  'z-z ) }  ♦ 

^  ( t: ) ( til) }  A/r/ 4  A//, ^ ( z,: )  V  r._, ) } 

*  -^zcz,.,)  {®4(^z-z)}-®^r  ’  ^Z(Z,.,)  {''('i^}'^Z(Zi_,)  {-^4  (*Z-z)}  A-ff^/ 

*^Z(.i,,){’’('z)'’^('z-z)}  rnf^y 
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Using  Eq  (3)  for  the  last  three  terms  we  get; 

The  middle  three  terms  are  easily  evaluated: 


{^*(h  )e*(Vi) 

}ab/  = 

(78) 

The  first  three  terms  require  an  expression  of  in  terms  of  £*(/; . We  will  use  a  similar 
development  to  what  we  did  in  Section  3.2.4  where  we  derived  an  expression  for  but  we  need 
6t(t;  )  for  this  expression.  We  start  with  the  definition  of  ej,t-)  and  decrement  one  time  sample: 

)  =  Xj-C  tj)  -  ^  ifC  ) 

=  I  , )  *  it  J.H  ( I ) .  Gj-w  ( /|._ ,  )  j  -  [  h- 1 )  * -B*  “  ( ] 

=  j.Xj,( /j.  I  )  -  (  4>  j.- A  ).ij(  <(.  ]  )  ♦  JBj.-  ABj^  j  j  u  ( *  Gj-w  (  j  ) 

=  1 )  -  ^* ( h’- 1 )  ]  ♦  4  *t^*(  </.  I ) ♦  ABjM  (/,.,)•  Gj.w  ( i,. ,  ) 

=  ®  _ ,  )  -  (  l,-.  1  )  -  i'jTj  (  ,  )  ]  ♦  A  /(.  1  )  *  JS'jrj(  l|.  ,  )  I 

»ABjj<(r,..,).Gj.H- (/._,) 

=  )]♦  A<l>t^j(r,'.,) 

♦  A  ®  j  z  (  _  j  )  -  B (  h  - 1  )  ]  *  ABj^ii  C  ^  }  f  G J.W  (_  ^  ) 

.  A  V  (  r,-_ ,)- (ffj.- ABt)i*(  <:.,)].  AB;t«  ( 1 ).  Gj-m- (/,.  1 ) 

~  ^  ® r'^*^r['*r(  ^i-i ) ” h-i  )] ■  ^  ^i  - 1  ^ h-i ) 

.  A  , (  q\ , ) .  A  4., IT* x^( (  r:. ,  )  1  ^  A  4., JS:* V  (  t, . ,  ) 

♦  A  4>j^j  ABj^j (  f(\  1 ) «  ABj^u  (  r, .  1 )  »  Gj.w  (  r._ , ) 

=  4.j.(/-j:jBj.)e*(r;_i)-(*r-A*k)^kV(/,.,)-(*j.-A4.t)r*ABt^*(r:.,) 

♦  A  (  'i  - 1  )  ^  A  ^  )  »  AB^ii  (  t,_  j  (  r,_  j  ) 
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(79) 


*  [  A  «■  t -  ]^* (  '/- 1 )  •  ( *, . , ) 

We  apply  this  equation  repeatedly  to  get: 

i-\ 


n-O 

/-I 


n  -  0 


(80) 


Note  that  we  have  separated  the  noise  terms  from  the  deterministic  terms  in  preparation  of  taking  the 
expected  value  of  this  equation. 

Now  to  use  this  in  the  first  term  of  Eq  (76): 

t-i 

*  BtJ:  [*r(^-^*^r)-A*t^*^r]"|(  a**- .„)£[(/, -.,)}|f/ 

B  -  0  ^  ■* 

/-! 

,  [<I.,(/-A:,^r).A4.,r*Fr]"|AB,£2(^  ,.„)£[(/, ■.,)}]ff/ 

B-  0  ‘  ' 

t-l 

.  ffrE  l<^r(^-J^kBr)^^<^,^tBrrjGrEz(i,  .)  { ’^  ( '.- 1 -b  )  }^Z(r,  ,){  e[(  »/-/)} 

n«0  ‘  ■' 

1-1 

n  -  0  ^  ^ 

=  .ff J.[ 4. J. ( / -  J.) .  A  4. , jr^ { e*( q. , ) e[(  r,:. ,)}b^^ 

i-i 

wO  ^  •' 

/-I 

♦  BjZ  [  *  A  j,j"f  ASjB  (/,._]  ejj,  (  r,_j)  I  \b-p 

-  Ar^[4>j.(/-z,fl'j.).A4.,Ar*ff^]'£^(^^){e*(t:,)eJ'(r:_,)}F/ 

1-1 

BtZ  [*r(^-'K'*^r)*A4.^i:,ffrr[(  a**- '^t^*A^)^*('i^i-B)^AB,«  (/..,.„)] 


n  =  0 
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Using  the  same  procedure  for  the  second  term  in  Eq  (76): 

£'z(^  (r,-.;)|A£ry 

i-i 

n  -  0 

•  )^z(^  j)|^t(^/-i  -  n)^t  (^1-/)}  I 

i-i 

n-O  ^  J 

/-I 

*  ^tY  ^  [^r^Z(^.,){  (  ^i-  1  -n)}^Z(,ti  i){^k  C  'i-i)}  1 

M-  0  ■  J 

Z-1 

n-O  ^  J 

^rf ^ *i^*^r]*^z(^.,){ *<-/)} ^*  (t.  ^)Affj. 
z- 1 

n-O 


Finally,  we  use  Eq  (79)  to  expand  on  the  third  term  in  Eq  (76) 

=  M*r(^-^*^r)-A-^.J*^*^rr^z«,,){^*(Vz)}^z(z,.){v^(^.z)} 

z- 1 

z- 1 

n-O  ^  J 

l-\ 

•BrT.  [t,(/-I.ffr)- i®.*'.ffrr|GA(, ,){>■(>, -.-,))«z(, ,){-’■«, .,))| 

n  —  0  ■* 

/-I 

•  BrT.  |*r('-*.ffr)-i®.r,ffrr(-®,i.£2(,.,,{v(l,_...)v'-(, ,.,))] 

n-O  J 

/-I 


60 


Note  that  6 («-/+!)  is  only  nonzero  for  «  =  M,  therefore  there  is  only  one  term  in  the  summation  of 


Eq  (84)  that  is  nonzero.  Also,  the  summation  is  only  valid  for  /- 1  > 0,  or  />  1 .  Thus : 


(84) 


We  put  Eq  (77),  Eq  (78),  Eq  (81),  Eq  (82),  and  Eq  (84)  into  Eq  (76),  to  get  an  expression 
for  the  residual  autocorrelation  with  mismodeling  in  the  output  matrix,  the  control  input  matrix,  and 
the  state  transition  matrix: 

l-\ 

n  -  0 

I- 1 

a  V*^rrl(  a**- (/,.,„)] 

n  ■  0 


This  is  a  ponderous  expression  to  evaluate  because  we  have  allowed  for  any  mismodeling  in 
the  control  input  matrix,  output  matrix,  and  state  transition  matrix.  Most  mismodeling  in  the  MMAE 
filter  banks  occurs  in  only  one  of  these  matrices.  Research  [11,  12,19, 41, 42, 43, 47, 48, 49,  52, 58, 
60,  61]  has  shown  that  the  MMAE  identifies  flight  control  sensor  failures  quite  well,  so  we  will 
target  our  research  to  improve  the  MMAE  performance  in  detecting  actuator  failures.  Thus,  we  will 
limit  the  mismodeling  for  the  rest  of  this  development  to  a  mismodeled  control  input  matrix.  With 
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this  assumption,  we  get  Aff  =  0  and  AO  =  0,  which  simplifies  the  previous  equation  for  the  residual 
autocorrelation  to: 


To  complete  this  derivation,  we  need  to  develop  an  expression  for  the  first  term  in  Eq  (86).  We  start 
with: 


=  I  ^  ^1-  / )  *  *i-  / )  ■  *1-  / )  ]  ^-z  ^ *1-1 )  ]  } 

'  { [*r(  *i-z )  ■  *i-z)  *  ^ ®*(  *i-z) ]  *1-1 )  ■  *1-1 )  ]  } 

°  ^z(^.i)  { [®r(*i-i)*  ^^*(*i-/)][^r( *<-/)*  ^^*(*i-i)]  } 

■  ^z(<i_,)  {  *1-1 )  *1-1 )  }* ■^z(^_,)  {  *i-i )  ^ *1-1 ) } 

^z(ti  ]){  ^^*(*i-i)®r(*i-i)}* ^z(ti  I )  {  ^ ( *1- 1 )  ^ ( *<■  -1 )  } 

=  *  ^Z(,ti_^)  {  ^r(*i-l)}  ^^*(*i-l)’  {  ^r(*i-l)} 

^^e4(<;_,)Aej(i;.,)*' 

=  Pj-  -t  A€^(ti_i)^ 


where  Eq  (27)  is  used  to  define  P^iti.;),  and  we  replace  this  with  the  steady  state  covariance  matrix 
P-f  in  the  last  line  of  Eq  (87).  We  can  use  this  expression  in  Eq  (86)  to  compute  an  estimate  of  the 
autocorrelation  of  the  residual  given  a  mismodeling  in  the  input  matrix,  and  then  use  Eq  (86)  to 
compute  the  mean  of  the  power  spectral  density  estimate  using: 
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where  we  can  use  the  relation  A/^il)  =  A^{-1)  since  the  residual  is  a  real  sequence.  We  can  use  Eq  (86) 
and  Eq  (87)  to  compute  the  residual  autocorrelation  matrix  and  put  this  into  Eq  (88)  to  get  the 
desired  mean  of  the  power  spectral  density  estimate. 

A  more  practical  method  of  computing  the  mean  of  the  power  spectral  density  estimate 
would  be  to  compute  the  mean  of  the  residual,  given  the  mismodeling  in  the  Kalman  filter,  using  the 
development  in  Section  3.2.4,  and  then  compute  the  power  spectral  density  estimate  using  Eq  (71)  or 
Eq  (72).  This  method  will  not  yield  the  correct  mean  of  the  spectral  density  estimate,  but  for  large 
residuals  it  will  provide  an  acceptable  approximation. 


develop  an  expression  for  the  covariance  of  the  power  spectral  density  estimate,  we  propose  a 
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vectorized  version  of  the  power  spectral  density  estimate  matrix  defined  in  Eq  (71).  To  form  this 
vector,  we  take  the  columns  of  the  ❖  matrix  and  stack  them  one  on  top  of  the  other.  Therefore: 


\  t,') 


’^3.1 


(89) 


where  m  is  the  number  of  elements  in  the  residual  vector  and  f  ]  is  the  ceiling  function,  which  rounds 
the  operand  to  the  nearest  integer  toward  positive  infinity.  Thus,  f  (/;  r,)  is  an  m^-dimensional 
vector  of  functions  of  frequency /.  This  construct  will  aid  in  notation  for  the  covariance  of  this 
estimate  since  the  covariance  will  be  a  matrix  of  scalars  using  this  vector  version,  versus  a  tensor  if 
we  use  the  matrix  version  of  the  power  spectral  density  estimate. 

Using  this  notation,  we  will  derive  an  expression  for  the  autocorrelation  of  the  power 
spectral  density  estimate.  We  will  derive  this  expression  for  a  single  element  of  the  autocorrelation 
matrix  using  a  and  b  to  denote  that  particular  element  (note  that  a  and  b  can  each  take  on  integer 
values  from  1  to  m^).  Also,  we  are  denoting  a  particular  element  of  the  residual  vector  using;,  k,  /,  or 
«  as  the  index. 


64 


N-\ 


E 

«.-(w-i) 


a,  ^ 

w-i-iji 


^  jr-0 


V-1 


W-l-|p| 


p=-(w-i) 


N 


c  =  0 


.  J/-1  W-1  AT-l-l^l 

—  E  E  E  E  ,){ '•-.('i-c-ii-i)} 

N  q^-{N-\)  p‘-{N-\)  j-0  c  =  0 

•  exp(-72  «/(/>«?)) 


where  j  =  —,  k  =  a  -  {j  -  \  )m  ,  I  =  —  ,n  =  fc-(/-l)f 

m  m 


(90) 


This  expression  involves  the  fourth  moment  of  a  Gaussian  random  variable,  which  is 
difficult  to  expand  if  we  allow  any  mismodeling  in  the  Kalman  filter.  We  intend  to  use  this  residual 
correlation  Kalman  filter  bank  with  the  standard  hypothesis  testing  algorithm  that  is  desaibed  in 
Section  3.3.1.  This  algorithm  only  requires  the  covariance  of  the  power  spectral  density  estimate  for 
the  case  of  no  mismodeling  in  the  Kalman  filter.  We  showed  earlier  that  this  produces  a  zero-mean 
residual  and  we  can  compute  the  steady  state  covariance  of  the  residual.  Therefore,  we  can  use  the 
standard  development  of  the  fourth  moment  of  a  zero-mean  Gaussian  random  variable  [e.g.  37:  107] 
to  expand  this  last  expression,  since  the  fourth  moment  can  be  expressed  in  terms  of  the  second 
moment,  which  is  the  covariance  of  the  residual.  Thus; 


1 


W-l 

-  y 

At  J--(W-I)  ;;.-(W-l)  x-O 


a,  b 

N_l  JV-Mg|  N-l-\p\ 


t-o 


■exp{-j2nf(p^q)) 


(91) 
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where  ,{s  - 1)  denotes  the  j,  I  th  element  of  theA^-C  s  - 1)  matrix. 

This  expression  is  very  computationally  intensive.  Note  that  we  are  computing  the 
covariance  of  the  power  spectral  density  estimate  for  the  case  where  there  is  no  mismodeling,  so 


Ar(t) 


f  A  j  (  0  ),  for  t  =  0 


0.  for  t  t  0 


(92) 


where  A^(0)  is  the  precomputed  steady  state  Kalman  filter  residua!  covariance  matrix  for  the  true 
Kalman  filter.  Thus,  we  can  approximate  Eq  (91)  by  checking  each  of  the  time  arguments  of  the 
residual  covariance  matrices  and  assume  that  the  residual  covariance  matrix  is  0  if  the  time  argument 
is  not  zero. 

Earlier  we  showed  that,  for  the  case  of  a  mismodeled  input  matrix,  corresponding  to  an 
actuator  failure,  the  residual  covariance  matrices  are  all  equal.  Thus,  the  covariance  of  the  spectral 
density  estimate  will  be  the  same  for  any  of  the  actuator  failure  hypotheses. 

3.3  Hypothesis  Testing  Algorithm 

In  Section  3.2  we  showed  that  the  Kalman  filter  residuals  within  an  MMAE  structure  are 
normally  distributed  random  vectors  with  known  (and  precomputable)  steady  state  covariances,  but 
different  means.  Therefore,  the  residuals  will  have  different  distributions  that  are  dependent  on  the 
hypotheses  of  their  internal  models.  The  various  hypotheses  that  were  explored  for  this  research  are 
defined  by  and  AO*.  We  will  denote  a  set  of  these  parameters  as  0  and  the  parameter 

space  of  all  possible  parameter  variations  as  9.  A  particular  hypothesis  is  constructed  by  defining  a 


subset  of  the  parameter  space: 


We  also  assume  that  the  various  hypotheses  form  a  disjoint  covering  of  the  parameter  space,  so  that 

e  -  00  u  e,  u  -  u  (94) 

If  we  test  hg  versus  hj  versus  ...  then  we  have  a  N-ary  hypothesis  test.  If  N=2  we  have  a  binary 
hypothesis  test.  The  primary  hypothesis,  which  is  the  assumed  true  hypothesis  at  the  time  of  the  test, 
will  be  denoted  hg.  The  other  hypotheses  will  be  called  the  alternate  hypotheses.  Thus,  a  binary 
hypothesis  test  is  a  test  of  the  primary  hypothesis,  hg,  versus  the  alternative  hypothesis,  h^.  N-ary 
hypothesis  testing  is  an  extension  of  binary  hypothesis  testing  by  simply  making  N-1  binary 
hypothesis  tests  between  the  primary  hypothesis,  hg,  and  the  N-1  alternate  hypotheses,  hj  through 
to  obtain  the  desired  N-ary  hypothesis  test.  This  is  sufficient  because  the  covering  in  Eq  (94)  is 
assumed  to  be  disjoint.  If  each  subspace  0*  contains  a  single  element  (i.e.  represents  a  single  flight 
control  failure),  then  the  hypothesis 

(95) 

is  called  a  simple  hypothesis.  If  it  contains  more  than  one  element,  then  it  is  called  a  composite 
hypothesis.  We  will  assume  that  the  hypotheses  we  are  researching  are  simple  hypotheses  (each 
subspace  contains  only  a  single  flight  control  failure  hypothesis)  and  denote  a  particular  hypothesis 
as: 

A*:  0  =  6*  (96) 

There  are  two  types  of  errors  that  can  occur  when  testing  hypotheses.  Type  I  errors,  called 
false  alarms,  occur  when  the  alternate  hypothesis  is  chosen  when  the  true  hypothesis  is  the  primary 
hypothesis.  The  probability  that  this  type  of  error  can  occur  is  called  either  the  probability  of  false 
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alarm,  denoted  Pp^,  or  the  size  of  the  test,  denoted  a.  Type  IT  errors,  called  missed  detection,  occur 
when  the  primary  hypothesis  is  chosen  when  the  true  hypothesis  is  the  alternate  hypothesis.  The 
probability  that  this  type  of  error  can  occur  is  called  the  probability  of  a  miss.  This  probability  is  1- 
Pp,  where  Pp,  is  the  probability  of  detection,  also  called  the  power  of  the  test. 

Please  note,  the  definitions  for  a  false  alarm  and  a  missed  detection  that  are  used  for 
Neyman-Pearson  hypothesis  testing.  Researchers  in  the  failure  identification  field  usually  define  a 
false  alarm  as  a  declaration  of  a  failure,  when  in  fact  one  had  not  occurred,  and  a  missed  detection  as 
a  declaration  of  no-failure,  when  in  fact  a  failure  had  occurred.  These  definitions  agree  with  the 
Neyman-Pearson  definitions  above,  if  the  current  hypothesis  is  the  no-failure  hypothesis. 
Disagreement  occurs  when  a  failure  hypothesis  is  in  force.  For  example,  if  the  assumed  primary 
hypothesis  is  a  failed  elevator,  and  the  hypothesis  test  declares  a  no-failure  when  in  fact  the  elevator 
failure  remains  true,  the  Neyman-Pearson  definition  would  define  this  a  false  alarm  but  the  usual 
terminology  would  define  this  as  a  missed  detection.  We  chiKise  to  use  the  Neyman-Pearson 
definitions  because  they  are  defined,  without  ambiguity,  by  the  design  parameters  Pp^  and  the  Pp, 

3.3. 1  Standard  Hypothesis  Testing  Algorithm.  First,  we  will  give  an  example 

to  demonstrate  the  conceptual  basis  of  the  Standard  Hypothesis  Testing  Algorithm  (SHTA),  and  then 
we'll  derive  the  specific  equations  of  this  algorithm.  Let  us  assign  Kalman  filter  1  the  fully  functional 
aircraft  model,  and  Kalman  filter  2  the  left  elevator  failure  model.  Each  of  these  filters  would  use  the 
commanded  input  and  its  internal  model  to  compute  an  estimate  of  what  the  measurements  should  be. 
Included  in  this  computation  is  a  statistical  model  of  both  the  noise  in  these  measurements  and  the 
disturbances  that  drive  the  system  dynamics.  The  Kalman  filter  algorithm  uses  these  models  to 
calculate  an  expected  value  and  covariance  matrix  of  the  measurement,  before  the  actual 
measurement  is  taken  at  a  particular  sample  time.  The  actual  measurement  is  then  used  to  form  the 
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residual,  which  is  the  difference  between  the  actual  measurement  and  the  computed  prediction  of  that 
measurement.  If  the  aircraft  has  no  failures,  then  the  residual  from  filter  1,  whose  internal  model 
hypothesizes  that  there  are  no  failures,  would  be  much  smaller,  relative  to  its  own  internally 
computed  covariance,  than  the  residual  in  filter  2. 

The  SHTA  would  use  the  information  in  the  residuals  from  the  filter  bank  to  compute  the 
relative  probabilities  of  each  of  the  hypothesized  models.  The  algorithm  would  assign  the  highest 
probability  to  filter  1  since  it  has  the  smallest  residual.  Now,  let  us  assume  that  a  left  elevator  failure 
occurs,  thus  the  residual  in  filter  2  would  become  quite  small  and  the  residual  in  filter  1  would  grow. 
The  hypothesis  testing  algorithm  would  tlien  assign  less  probability  to  the  hypothesis  of  filter  1  and 
more  to  the  hypothesis  of  filter  2. 

3. 3. 1.1  Basic  Equations.  The  Standard  Hypothesis  Testing  Algorithm 

(SHTA)  simultaneously  tests  the  residuals  of  the  Kalman  filter  bank  under  multiple  hypotheses.  In 
Section  3.2.4. 1  we  found  that  if  the  Kalman  filter  model  matches  the  true  system  model,  the  residual 
has  a  mean  of  zero,  and  in  Section  3.2.3  wc  found  that  the  residual  will  have  a  precomputable 
covariance  matrix,  A^.  Therefore,  this  Kalman  filter  residual  is  a  white  Gaussian  sequence  of  mean 
zero  and  covariance 

\  *  R,  (97) 

Since  we  know  that  the  residual  is  a  Gaussian  vector  with  zero  mean  and  covariance  A;^,  we  substitute 
these  values  into  the  known  expression  for  a  Gaussian  conditional  density  function.  Therefore  we  get 
that  the  conditional  density  function  of  the  measurement  ( z  )  at  t,  for  the  kth  Kalman  filter, 
conditioned  on  the  measurement  history  (  Z(/,  . ,)  =  [  z^(r,  )•  ■••■  z^(t,  .  ,)]^ ),  is 
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where  P  j  = 


(98) 


^2(,t,)\h^.z(t^_^')  I  "  P^expl-} 

- -i - - 

(2  k  Mil  '  ’ 


The  scalar  likelihood  quotient  will  be  defined  as; 

=  r/’(t,)^t’rj(q)  (99) 

If  the  RCKFB,  developed  in  Section  3.2.6,  is  implemented  with  the  Standard  Hypothesis 
Testing  Algorithm  (SHTA),  we  modify  Eq  (98)  and  Eq  (99)  to  use  the  power  spectral  density 
estimate  from  the  RCKFB.  Two  versions  of  the  power  spectral  density  estimate  were  developed  in 
Section  3.2.6,  each  of  which  requires  a  slightly  different  modification  to  Eqs  (98)  and  (99).  The 
periodogram  version,  Eq  (7 1 ),  produces  a  matrix  of  power  spectral  density  estimates  because  it 
included  the  aoss  power  spectral  density  estimates  between  different  elements  of  the  residual.  We 
computed  the  mean  of  this  estimate,  Eq  (88)  with  4^(0  given  by  Eq  (85),  and  then  formed  a  vector 
version  of  this  matrix,  Eq  (89),  to  help  in  developing  an  expression  for  the  covariance  of  the 
estimate,  Eq  (91).  We  introduced  a  simpler,  less  accurate  version  of  the  power  spectral  density 
estimate  in  Eq  (72).  This  version  produces  a  vector  because  it  does  not  compute  the  cross  power 
spectral  density  estimate  between  the  elements  of  the  residual,  so  it  provides  a  subset  of  the  estimates 
that  are  computed  using  Eq  (71).  The  mean  for  this  Fourier-transform-based  version,  can  be 
computed  either  by  using  the  appropriate  terms  of  Eq  (88),  or  by  using  the  computation  of  the  mean 
of  the  residual,  developed  in  Section  3.2,4,  in  place  of  the  residual  in  Eq  (72).  The  covariance  of  the 
estimate  can  be  computed  using  the  appropriate  terms  of  the  covariance  matrix  from  Eq  (91). 
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Another  method  would  be  to  take  several  samples  of  the  power  spectral  density  estimate,  subtract  the 
mean  from  the  samples,  and  estimate  the  covariance  of  the  samples.  Many  software  packages,  such 
as  MATLAB,  have  an  intrinsic  sample  covariance  computation  available,  which  makes  this  method 
quite  easy  to  implement. 

In  practice,  the  mean  for  the  no-failure  hypothesis  is  very  small  when  compared  to  the  mean 
of  the  spectral  density  estimate  for  an  actuator  failure  hypothesis.  This  is  clearly  seen  in  Figure  8, 
where  the  no-failure  case  has  a  mean  of  about  1,  while  the  failed  case  has  a  mean  of  about  40.  Since 
the  disparity  between  the  hypotheses  is  so  great,  we  choose  to  mode)  the  power  spectral  density 
estimate  as  a  zero-mean  process  with  a  covariance  matrix  that  is  assumed  to  be  constant  and 
computed  using  400  samples  of  the  fully  functional  residual  power  spectral  density  estimate.  We 
evaluate  t  (/;  q)  at/ -f^,  where /q  is  the  known  frequency  of  the  input.  This  gives  us  the  modified 
version  of  Eq  (98): 


=  P^expl-l 

y^here  p  =  - -1 - —  ^  1  •  1  =  (  -  ^  ;  /  )4  (4  ; \ 


(100) 


and  Eq  (99)  becomes; 

9^(/)=  #*(/o;/)>4^Vt(/o;/) 


(101) 


At  this  point,  the  development  of  the  SHTA  for  the  Standard  Kalman  Filter  Bank  (SKFB),  Eqs  (98) 
and  (99),  parallels  the  development  of  the  SHTA  for  the  RCKFB,  Eqs  (100)  and  (101).  We  present 
modifications  to  the  SHTA  for  the  SKFB,  with  the  understanding  that  the  same  development  would 
be  used  for  the  SHTA  for  the  RCKFB,  including  the  "P  Stripping"  modification  discussed  in  Section 
3. 3. 1.2,  which  may  aeate  numerical  conditioning  problems  due  to  the  high  dimensionality  of 
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We  can  use  these  equations  as  a  relative  measure  of  the  accuracy  of  the  Kalman  filter  model. 
In  Section  3.2.4  we  found  that  if  the  Kalman  filter  model  has  some  mismodeling  errors,  the  mean  of 
the  residual  will  be  nonzero,  and  if  the  Kalman  filter  model  is  accurate,  the  mean  of  the  residual  is 
zero.  The  Kalman  filter  that  is  most  accurate  will  have  a  scalar  likelihood  quotient  (the  norm  of  the 
residual,  scaled  by  tlie  inverse  of  the  precomputed  covariance  matrix)  that  is  smaller  that  the  scalar 
likelihood  quotient  from  a  Kalman  filter  with  mismodeling  errors.  Thus,  Eq  (99)  becomes  a 
quadratic  penally  for  any  mismodeling  in  the  Kalman  filter  model.  Conversely,  if  we  use  the  scalar 
likelihood  quotient  as  a  negative  exponent,  as  in  the  "dot"  term  of  Eq  (98),  we  will  have  an 
exponential  indicator  of  the  accuracy  of  the  Kalman  filter  model.  Note  that  if  the  residual  is  larger 
than  expected  due  to  modeling  errors,  the  scalar  likelihood  quotient  is  large  and  the  "dot"  term  will  be 
small.  If  the  "dot"  term  is  scaled  by  the  P*  term,  we  have  a  conditional  probability  density. 

We  define  the  conditional  probability  for  a  particular  hypothesis  as: 

/>*(',)  =  prob  I  />  =  Aj  I  Z(t,)  =  z.  }  (102) 


We  can  compute  the  conditional  probability  for  a  particular  hypothesis  by  comparing  the  conditional 
probability  density  for  the  current  measurement,  assuming  that  particular  hypothesis,  to  the  densities 
associated  with  all  the  other  hypotheses.  The  conditional  probability  of  a  specific  hypothesis  can  be 
shown  [34]  to  be: 


^  ^  fz(t,)\h.Z(,t,  .)  i  ^-1  ) 

Pk(t^)  =  — - : - 

E  /i(r,)lA.2((,_,)  (2.1  V •  Ty(f.-i) 


(103) 


In  this  equation  we  use  the  prior  conditional  probabilities,  Pj.(t, . ;),  to  weigh  the  conditional 
probabilities  of  the  current  measurements,  assuming  each  hypothesis,  and  tlien  normalize  it  over  the 
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complete  set  of  such  numerator  terms.  For  failure  identification  applications,  where  we  usually  want 
to  choose  the  most  likely  hypothesis  out  of  the  set  of  possible  hypotheses,  we  would  choose  the 
hypothesis  with  the  largest  conditional  probability.  These  conditional  probabilities  can  also  be  used 
to  weight  and  blend  the  various  hypotheses,  depending  on  the  particular  application. 

In  practice,  these  conditional  probabilities  will  fluctuate  rapidly  from  one  time  sample  to  the 
next.  If  the  hypothesis  testing  algorithm  is  choosing  the  hypothesis  with  the  highest  conditional 
probability,  these  fluctuations  could  cause  momentary  incorrect  hypothesis  declarations.  To  alleviate 
this  phenomenon,  the  conditional  probability  calculations  are  modified  (as  described  in  the  next  few 
sections)  and  compared  to  a  decision  threshold.  The  decision  threshold  is  defined  such  that  a 
hypothesis  is  chosen  only  if  its  conditional  probability  is  greater  than  the  threshold.  The  threshold  is 
set  to  avoid  momentary  incorrect  hypothesis  declarations  while  hopefully  providing  adequate 
hypothesis  testing  performance.  The  performance  is  usually  measures  by  the  lime  it  lakes  to  detect 
and  declare  tlie  correct  flight  control  failure  for  the  failure  identification  application.  Various 
modifications  to  the  conditional  probabilities  are  explained  below. 

3. 3. 1.2  P  Stripping.  Stevens  (43,  60]  found  that  certain  performance 

problems  could  be  reduced  by  modifying  Eq  (98).  He  altered  the  conditional  density  function  in  Eq 
(98)  by  removing  the  term,  which  was  used  to  make  the  area  under  the  density  function  equal  to  1. 
Thus  Eq  (98)  becomes: 

,)  (Zf)  = 

where  { •  }  =  j  -  i  r*  ( /,  )  ^^  '•*(0)} 

whereP^^  (z,.)is  the  unnormalized  conditional  probability  of  hypothesis Eq(103) 
becomes: 
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Pk(t.)  = 


(Zi)  •  Pk^h-\') 


(105) 


52  ^ hj.zit,  ’  Pj^U-\^ 

J-1 

If  all  tlie  scalar  likelihood  quotients  were  approximately  the  same  size  for  all  elemental 
filters,  Eq  (98)  and  Eq  (103)  would  put  the  highest  probability  on  the  elemental  filters  with  the 
smallest  jA^I  value.  This  is  an  inappropriate  weighting  since  the  size  of  lA^j  has  nothing  to  do  with 
the  correcmess  of  the  hypothesis  in  matching  the  current  real-world  failure  status.  Since  sensor 
failures  exhibit  themselves  as  a  row  of  H  going  to  zero,  filters  based  on  the  hypothesis  of  a  failed 
sensor  tend  to  have  smaller  lA^I  values,  and  thus  the  MMAE  will  be  prone  to  false  alarms  on  sensor 
failures. 

Eq  (104)  and  Eq  (105)  function  properly  with  the  P*  term  removed  because  the  denominator 
in  Eq  (105)  is  the  sum  of  all  possible  numerators,  so  the  probabilities  (p*)  will  still  sum  to  1  even  if 
the  area  under  each  of  the  modified  "densities"  of  Eq  ( 104),  which  are  the  densities  of  Eq  (98)  with 
Pi  removed,  is  no  longer  unity. 

3.3.1 .3  Lower  Bounding  of  the  Conditional  Probabilitie.s. 

Previous  research  [42,  43,  52,  60,  6 1  ]  also  found  that  a  lower  bound  needed  to  be  placed  on  the 
probabilities.  The  purpose  of  the  MMAE  is  to  make  quick  and  accurate  failure  identifications,  and  it 
was  found  that  if  some  of  the  probabilities  were  allowed  to  gel  loo  small,  it  took  a  long  time  for  the 
probability  in  the  correct  filter  to  grow  when  a  change  in  the  failure  status  occurred.  This  was  due  to 
the  fact  that  the  previous  conditional  probability,  Pt(t,.,),  was  so  small  for  the  new  correct  filter  model 
and  so  large  for  the  old  correct  model,  that  Eq  (105)  had  to  be  iterated  several  times  before  the  values 
would  change  significantly. 
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This  method  was  developed  as  an  ad  hoc  was  to  allow  for  time  varying  model  parameters.  If 
we  assumed  that  the  model  parameters  were  not  constant,  we  would  need  to  set  up  a  filter  bank  in 
which  each  filter  was  based  on  a  particular  time  history  of  the  model  parameters.  This  causes  a 
geometric  growth  of  the  filter  bank.  A  number  of  researchers  have  pursued  such  MMAE  designs,  as 
particularly  with  Markov  process  models  for  probability  variations  in  time  [13,  14,  32,  33],  but  some 
means  of  pruning  or  merging  branches  in  the  resulting  decision  tree  structure  is  needed  to  prevent  the 
impractical  growth  of  the  filter  hank.  Lower  bounding  the  conditional  probabilities  in  Fq  (105) 
allows  a  change  in  model  parameters  to  be  reflected  in  the  conditional  probabilities  without  such 
dramatic  bank  growth.  Menke  [47, 48, 49]  and  Stratton  [61]  found  that  0.001  was  a  good  lower 
bound  on  the  conditional  probabilities  for  an  application  using  about  ten  elemental  filters. 

3. 3. 1.4  Filter  Tuning.  Previous  experiments  [  19, 41]  found  that  an  incorrect 

hypothesis  may  be  momentarily  chosen  because  the  scalar  likelihood  quotients  representing  the 
incorrect  hypothesis  may  temporarily  be  smaller  than  the  likelihood  quotients  from  the  TOrrect 
hypothesis.  One  cause  of  this  phenomenon  was  the  norm  of  A,,  for  these  incorrect  hypotheses.  These 
hypotheses  had  smaller  values  for  the  norm  of  4*,  and  thus  had  smaller  scalar  likelihood  quotients  for 
a  given  vector  residual.  As  noted  in  Section  3.3. 1.2,  the  norm  of  4*  has  nothing  to  do  with  the 
correctness  of  the  hypothesis,  so  these  small  norm  values  were  causing  momentary  incorrect 
hypothesis  declarations. 

This  effect  was  corrected  by  tuning  the  Kalman  filters  that  used  the  incorrect  hypothesis 
models.  The  tuning  was  accomplished  by  increasing  the  modeled  value  of  measurement  noise,  R, 
which  is  called  adding  measurement  pseudonoise.  Adding  tlie  measurement  pseudonoise  changes  the 
Kalman  filter  gain  so  that  when  the  state  estimates  are  updated,  the  measurements  are  not  weighted 
as  heavily  as  they  were  before  the  pseudonoise  was  added.  This  technique  is  especially  appropriate  if 
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one  of  the  scalar  measurements  is  particularly  noisy  with  respect  to  the  other  scalar  measurements. 
Another  method  of  increasing  the  norm  of  is  to  increase  the  modeled  value  of  the  dynamic  noise 
covariance,  Q.  This  also  changes  the  Kalman  filter  gain,  but  with  this  technique  the  measurements 
are  weighted  more  heavily  then  the  predictions  of  those  measurements  based  on  the  Kalman  filter 
model.  If  the  filters  are  mistuned,  the  accuracy  of  the  state  estimates  is  much  poorer  than  can  be 
achieved  with  better  tuning  of  the  Kalman  filter  gain  values.  The  appropriate  amount  of  pseudonoise 
had  to  be  found  through  extensive  simulation  to  prevent  performance  degradation  due  to  mistuning. 

3.3. 1 .5  Smoothing  of  the  Conditional  Probabilitie.s.  The  false  alarm 

rate  in  some  previous  research  [42, 43, 52, 60, 61]  was  found  to  be  too  high  in  most  simulation 
examples.  To  decrease  the  false  alarm  rate,  the  conditional  probabilities  were  averaged  over  several 
time  samples.  A  set  of  N  conditional  probabilities  was  collected  and  then  averaged  using: 


n-\ 


^  E  Pjih-j) 


N 


j~o 


(106) 


This  average  was  then  compared  to  the  threshold  and  a  hypothesis  was  chosen.  This  decreased  the 
effects  of  momentary  false  alarms,  but  it  also  delayed  the  hypothesis  identification. 

3.3. 1.6  Exponential  Penalty  Increase.  The  SHTA  was  found  to  be  much 

slower  in  declaring  the  correct  failure  for  certain  failure  hypotheses  than  for  others.  These  failures 
took  much  longer  to  identify  because  the  effects  of  the  failure  had  to  grow  for  some  time  before  they 
were  disambiguated  from  other  failure  hypotheses.  To  remedy  this,  a  scalar  multiplier  was 
substituted  into  Eq  (104)  to  get: 

<•!  =  (  -  r  r/(/,)4t  Vj(/.)  }  (107) 
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Research  [19, 41]  has  shown  that  by  increasing  the  exponential  penalty  for  an  incorrect 
hypothesis,  that  is  increasing  F  in  Eq  (107),  caused  the  SHTA  to  choose  the  correct  hypothesis 
faster.  However,  the  number  of  false  alarms  also  increased  when  this  penalty  was  increased.  The 
penalty  needed  to  produce  the  desired  performance  had  to  be  found  through  extensive  simulation. 

3.3. 1.7  Propagating  Several  Time  Samples  before  Updating. 

Another  approach  to  rectifying  the  problem  presented  in  the  last  section  is  to  propagate  the  state 
estimate  several  time  periods  before  updating.  Figure  10  shows  graphically  how  this  approach  would 
work  with  a  scalar  state  estimate.  Normally,  the  state  estimates  are  updated  at  each  time  sample, 
which  moves  the  state  estimate  closer  to  the  measurement,  as  shown  in  Figure  10.  If  the  state 
estimate  is  propagated  without  update  for  several  time  samples,  the  distance  between  the  state 
estimate  and  the  measurement  will  be  much  greater,  if  the  hypothesized  model  used  by  the  filter  is 
incorrect.  This  is  clear  in  Figure  10  where  thei:,,  state  estimate  is  propagated  without  update  for 
three  iterations  and  the  distance  between  ^ ;e,(//,j)andz(/,.,3)  is  much  larger  than  the  distance 
between  ( q,, )  and  z(t,+3).  This  will  produce  a  larger  residual  ,  which  produces  a  larger  penalty 
for  mismodeling,  thus  helping  to  disambiguate  the  hypothc.scs.  Others  [20]  have  researched  this 
approach. 

3.3.2  Neyman-Pearson  Based  Hypothesis  Testing  Algorithm. 

3.3.2. 1  Basic  Equations,  in  this  section  we  will  develop  an  alternative 

hypothesis  testing  algorithm  based  on  the  Neyman-Pearson  Lemma.  We  will  first  define  a  binary 
hypothesis  test  and  then  present  the  Neyman-Pearson  Lemma.  Then  we'll  show  that  this  alternative 
hypothesis  testing  algorithm  requires  only  a  single  Kalman  filter  residual,  instead  of  a  complete  bank 
of  Kalman  filters. 
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Figure  10  -  Comparison  of  Updated  Slate  Estimates  with  Nonupdated  State  Estimates. 


A  binary  test  of  hg.  0  e  Qq  (the  parameter  sol  is  in  the  parameter  subspacc  0„)  versus  h,\  9 
e  0,  of  a  residual  would  take  the  form 


1 1  -  A,,  r  6  R 

I  0  -  Ap,  r  e  A  . 


(108) 


This  is  read  "the  test  function  (t)(r)  equals  1 ,  and  hypothesis  /t,  is  accepted,  if  the  residual  lies  in  the 
hg  rejection  region  R.  The  test  function  equals  zero,  and  hypothesis  hg  is  accepted,  if  the  residual  lies 
in  the  hg  acceptance  region  A." 

Considering  the  flight  control  failure  detection  application,  our  desire  is  to  define  the 
acceptance  and  rejection  regions  such  that  the  detection  probability,  is  maximized  for  a  given 
probability  of  false  alarm,  Pp^.  We  can  usually  accept  a  certain  rate  of  false  alarms,  but,  it  is 
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essential  not  to  miss  detecting  the  failure  with  this  particular  application.  There  are,  however,  other 
applications  that  might  need  to  minimize  the  probability  of  a  false  alarm  for  a  given  detection 
probability. 

Now  that  we  have  defined  the  basic  binary  hypothesis  test,  we  can  introduce  the  Neyman- 
Pearson  Lemma  [53;  107]: 

Nevman-Pearson  Lemma 

Let  0={  0n,  0,  }  and  assume  that  r  has  a  density  or  probability  mass  function,  denoted /g(r). 
The  subscript  denotes  that  the  function  exists  over  the  whole  parameter  space  0.  This  function  will 
be  subscripted  with  0q  or  0)  to  denote  the  function  for  a  particular  subset  of  tlie  parameter  space. 

The  test  of  the  form 


1.  '//«('■)>  k/gtr) 


r  }  =  j  if fQ^ir)  = 

0.  i//6.(0  <  k/e  (O 


(109) 


for  some  test  threshold  k  0  and  some  0  <  C  ^  L  is  the  most  powerful  test  (die  test  with  the  largest 
Pq)  of  size  a  >  0  (for  a  given  Pp^)  for  testing  hg-.  0  =  0o  versus  li,:  0  =  0,.  These  tests  are  unique 
except  perhaps  on  sets  of  probability  zero  under  /(„  and  /;,.  When  (})  (r)  =  1  we  choose  hj,  and  when 
(|)  (r)  =  0  we  choose  hg.  When  cfi  (r)  =  C  we  select  /i,  with  probability  C- 

This  lemma  can  be  rewritten  by  defining  the  likelihood  ratio  as; 

/(r)  =  - 

(110) 
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so  that  the  most  powerful  test  of  size  a  for  testing  hg.  8  =  0o  versus  h,:  0  =  0j  is  a  likelihood  ratio 
test  of  the  form 


1.  i//(r)>k 

4>  {»•)  =  -  C.  i/l(r)  =  k 

0,  i/l(r)<k  (111) 

If  I  {r)-k  with  probability  zero,  then  C  =  0  and  the  threshold  k,  for  a  given  size  a ,  can  be  computed 
from: 

“  =  -Pej  1  =  l\(l)dl 

t  (112) 

where  ( O  is  tlie  density  function  of  /  (r)  under  hg.  The  power  of  the  test,  P^,  can  be  obtained 
from: 

^  -  Pd  =  [fa,indl. 

k  (113) 

We  will  now  show  that  the  likelih(xxl  ratio  in  Eq  (1 10),  which  operates  on  a  single  residual, 
is  equivalent  to  a  likelihood  ratio  that  operates  on  two  different  residuals.  First,  we  assume  that  the 
density  functions  used  in  the  likelihood  ratio  are  conditioned  on  the  measurement  history  Z(r,_,).  We 
will  compare  two  residuals  from  two  different  Kalman  filters,  rj  and  r*.  The  Kalman  filter  models 
used  in  these  filters  use  two  different  hypotheses  that  differ  only  in  the  control  input  matrix,  output 
matrix,  or  dynamics  matrix.  We  will  compare  tlie  two  likelihood  ratios  to  check  for  equality: 

/e/'-*)  '  /e/'-p 
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(114) 


Po  {  -^['’k-^Z(t,_^){'’k}\h^f^0  '['■*-'^Z(^.l){''*}l/.„]} 

7  Pi  exp  j  -^[0-£z(<i..)h}lsf-^i*[0-^z(<,..){0}ll..]} 


For  simplicity  of  notation,  we  are  suppressing  the  time  argument,  t,,  of  the  residuals.  Examining  Eq 
(1 14),  we  note  that  the  numerators  and  leading  coefficients  are  equal.  Thus  the  equality  will  hold  if 
and  only  if: 


'■*  -  -^ZCl,., ){'■*}!*„  =  '■j  -  ^Z(l,.,)('‘j}lA, 


(115) 


We  assume  that  the  hypothesis,  /if,,  can  consist  of  any  of  the  modeling  differences  that  we  are 
investigating,  namely  AB,  AH,  and  A^>.  We  substitute  Eq  (35)  and  Eq  (36)  into  Eq  (1 15)  and 
cancel  terms  to  get: 


^r^r[  ^*(  h  - 1 )  ■  ^z(i,., ){ ( *1- 1 ) }  ]  '  \ )  *  '’r^  ^ ) 

=  ®  rl  ^;(  h  - 1 )  ■  ^z(t,., )  {  ^;  (  h  - 1 ) }  |  *  H  j.  w  t.  ■^)  *  Vj,(  /^ ) 


”*  1  )  ■  ^2(1,., )  I  (  ^- 1  )  }  =  ^y(h-l)  '^Z(^,,  )  {  (  ^1- 1  )  } 


(116) 


We  then  expand  tliese  terms  using  the  definition  from  Eq  (23)  to  get: 

[  •*r(^-l  )  ]  ■  ^Z(i;,j){-’^r(^-l  )‘^t(h-i )} 

=  [  *r(^'i  - 1  )■  - 1 )  ]  ■  ■^z(t,.,){*r( *1-1 )" ■^j(  ii-i  ) } 

•*r(*i_i )  -  ■^fc(^i-i )  -  ^z(r,,,){'*7'(*i-i ’  ■^*(*1-1 ) 

=  *r(*i-i)  ■  ■^y(*i-i)  ■  ^z(i,.,){*r(*i-i )}  ’  ^y(*i-i) 

•*r(  *1-1 )  -  ^z(tj.,){-’*r(  *1- 1 *1-1 )  ■  ^Z(rj_i){*r(  *1-1^} 
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Therefore 


/a/'"*)  ■4o^''P  (118) 

Thus  we  can  perform  the  hypothesis  test  using  a  single  Kalman  filter  residual,  testing  it  for  multiple 
hypotheses,  and  this  is  equivalent  to  performing  the  hypothesis  test  using  multiple  residuals,  each 
one  tested  for  a  single  hypothesis. 

3.3. 2. 2  vSingle  Time  Sninple  Hypothesis  Testing.  Now  we  can  apply 

the  Neyman-Pearson  hypothesis  test  to  the  residual  of  a  Kalman  filter.  We  need  an  estimate  of  the 
residual  covariance  for  a  particular  hypothesis,  which  we  can  precompute,  as  described  in  Section 
3.2.1,  by  using  a  Kalman  filter  that  uses  the  hypothesized  system  model.  We  will  use  the  steady  state 
Kalman  filter  estimate  of  the  residual  covariance  matrix,  thus  the  residual  covariance  will  be 
considered  constant.  Note  that  the  estimate  for  the  residual  covariance  is  not  a  function  of  the  system 
input  matrix,  B.  Thus,  if  we  only  have  modeling  differences  in  the  system  input  matrix,  so  only  AB 
is  nonzero,  then  the  residual  covariance  matrix  will  be  identical  for  each  of  the  hypotheses.  We  will 
develop  the  Neyman-Pearson  based  hypothesis  testing  algorithm  for  the  case  of  AB  *  0,  which  is  an 
actuator  failure  for  the  flight  control  application.  Later  we  will  look  at  the  case  where  the  other 
modeling  differences  are  nonzero. 

Under  these  assumptions,  we  have  a  Neyman-Pearson  detector  for  "common  covariances, 
uncommon  means"  as  developed  by  Scharf  [53;  111].  We  will  develop  the  Neyman-Pearson  detector 
in  the  same  maimer  as  Scharf  does,  but  we  will  first  extend  it  to  a  single-time  sample  hypothesis  test 
using  a  multidimensional  random  variable  (a  single  time  sample  of  the  residual).  In  Section  3. 3.2. 3, 
we  will  extend  the  development  to  use  several  time  samples  of  the  residual. 
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We  will  start  with  a  single  time  sample  of  a  Kalman  filter  residual  and  design  a  binary  test 


using  the  Neyman-Pearson  Lemma.  We  first  need  to  compute  the  likelihood  ratio.  We've  shown  that 
the  residual  is  normally  distributed  with  covariance irrespective  of  the  hypothesis  that  is  currently 
in  force  in  the  detection  algorithm.  These  residual  covariances,  which  we  will  denote  for  the 
various  hypotheses  are  identical  only  if  the  modeling  differences  are  limited  to  AB  *  0  and  all  other 
differences  (i.e.,  AO  and  A//)  are  zero.  We  will  denote  the  mean  of  the  residual  given  the  past 
measurement  history  and  a  pariicular  hypothesis  as 

I*.  =  (119) 


Using  Eq  (1 10),  the  likelihood  ratio  is  a  ratio  of  two  normal  density  functions,  which  becomes  the 
ratio  or  two  exponential  tcrin.s  when  the  leading  P  terms  are  canceled.  Thus  the  likelihood  ratio 
becomes: 


/(/■(/.))  =  exp  |-  j[r(r,)-m,(/,)]^^*''[r(/,.)-mj(r,.)] 


=  exp  y  I  r  ( r,. )'  ‘  r  ( r. )  -  r  ( /. )'  A^-‘  iWj  ( r, )  -  (  q  r  ( r, )  *  iWj  ( t. )'  A^  ‘  (  q. ) 

=  [  »•  ( ti  ‘  "*1  ( t, )  ♦  Wi  (  <,  t  ’  r  (  q  )  -  wtj  (  t,-  Wj  ( t J 


2 

^r  .  -1 


sr  -1 


T  A  -I 


Note  that  each  term  within  the  exponential  is  a  scalar  and  we  can  use  one  of  the  properties  of  the 
covariance  matrix  to  get: 
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Y  =  m^it^  {A  ^  f  r  {t^)  =  '  r  ( /,. ) 

=  (»■(«,• /^*'' Wo  (',)  )^=  =  Wo(',)^^*' '•('/)  (121) 

Also,  we  can  use  the  natural  log  of  the  likelihood  ratio  in  place  of  the  likelihood  ratio  because  the 
natural  logarithm  is  a  monotonic  function.  Therefore  we  get  the  log  likelihood  ratio; 

L(r(/.))  =  ln{/(r(/,.))} 

=  Wi  (  tifA^  t^A^  ^  ('"i^  '"i(  ‘i  )  '  ”(,(  '”o(  )  ) 

=  ( Wi(  ti )  -  Wo(  ti )  Y^Y  •'(  'i )  -  ^  ( '”i(  )  -  "*o(  h  )  *■’  ( "‘iC ', )  *  Wo(  /, ) ) 

where  s^t^)  =  ( WjC /, ) -  w,(  t,  )  r ( /, ) 

and  b(/.)  = 


We  call  s(ti)  the  signal  and  b(r,)  the  bias  of  the  log  likelihood  function.  One  structure  for 
implementing  Eq  (122)  is  shown  in  Figure  11. 

Rewriting  the  Neyman-Pearson  Lemma  to  use  the  log  likelihood  function  in  Eq  (122),  we 
get: 


4- {»■(/,.)  } 


E  •^'('■(^i))  >11  =  hi(k) 
0,  L  (  r  (  )  )  s  n  . 


(123) 
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Figure  11  -  Implementation  of  the  Neyman-Pearson  Based  Hypothesis  Test 

Now  we  only  need  to  compute  the  threshold  ri  in  order  to  define  the  Neyman-Pearson 
hypothesis  test  algorithm  completely.  We  first  need  to  find  the  probability  density  function  for  the 
test  statistic  L{  rf/,. )).  We  note  the  L  is  a  linear  transformation  of  r,  therefore  L  is  also  normally 
distributed  and  we  need  only  compute  the  mean  and  covariance  to  define  the  probability  density 
function.  We'll  start  by  computing  the  mean  under  both  the  primary  and  alternate  hypotheses. 

First  we  define: 

so  that  the  test  statistic  becomes: 

=  ■»(<.) -h(tj).  ^24) 

Now  we  compute  the  mean  of  the  test  statistic,  given  that  the  primary  hypothesis  is  in  force: 
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=  dit,  f  A  *->  m„(  t^)-jdi  t.  fA  j  1  m,(  /. )  -  i  J  (  t,  f  A  *  '  m„  (  / . ) 

=  I  '^(f,)'’^t‘'('«o(^)-«i(^))  =  -T =  -7^(',) 

2  <*Vo<  i,/2  2  (125) 


where  D{t. )  will  be  called  the  hypothesis  discrimination  measure. 

Likewise,  we  compute  the  mean  of  the  test  statistic  given  that  the  alternative  hypothesis  is  in 

force: 

-  d  it.  fA,'  m,i  /, )  -  i  (  /,  fA,'  «.j(  r, )  -  i  rf  (  r.  (  /^  ) 

-\dit^fA,'[m,it^)-m,it,))-\ditjA,'dit,)-\Dit,). 

2  '  '2  2  1126) 


Now  to  compute  the  conditional  covariance  of  the  test  statistic  under  either  hypothesis: 

-  [^z(t,.){^(^)}^-2^za,,){-(h)}b(^)-b(r.)^] 

Now  use  Eq  (124)  to  get: 
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=  dit^fA-^d{t^)  =  D(/,). 


<^(/,) 


(128) 


We  assumed  that  the  residual  correlation  matrix  was  identical  for  each  of  the  hypotheses 
when  we  formed  the  likelihood  ratio  in  Eq  (120).  This  occurs  when  we  only  allow  differences  in  the 
control  input  matrix  (B),  which  corresponds  to  actuator  failures  for  flight  control  failure 
identification  applications.  We  noted  earlier,  that  model  differences  in  the  output  matrix  {H)  and 
state  transition  matrix  (O)  will  produce  different  residual  correlation  matrices.  We  will  now  look  at 
the  case  where  we  are  modeling  a  flight  control  sensor  failure. 

A  flight  control  sensor  failure  is  modeled  by  zeroing  a  row  of  the  output  matrix  that 
corresponds  to  the  particular  failed  sensor.  Thus,  the  element  of  the  vector  H  x{t;)  that  corresponds 
to  the  failed  sensor  will  be  zero.  Therefore,  this  element  of  the  residual  vector  will  primarily  affect 
the  row  and  column  of  the  residual  covariance  matrix  that  corresponds  to  the  failed  sensor  and  the 
other  elements  of  the  covariance  matrix  may  be  slightly  affected.  It  has  been  observed  that  the  row 
and  column  of  the  residual  covariance  matrix  that  correspond  to  the  failed  sensor  increase 
significantly,  which  is  appropriate  since  this  indicates  that  the  Kalman  filter  algorithm  based  on  such 
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a  residual  covariance  matrix  will  not  rely  on  the  failed  sensor  measurement  as  much  as  it  did  before 
the  failure.  One  approach  to  detecting  flight  control  sensor  failures  would  be  to  conduct  a  Neyman- 
Pearson  hypothesis  test,  assuming  that  the  residual  correlation  matrices  are  equal  to  the  computed 
residual  correlation  matrix  for  the  sensor  failure  (and  correspondingly  for  the  element  of  H  x{t;) 
corresponding  to  the  failed  sensor  to  be  set  to  zero)  This  approach  will  work  well  if  the  effect  of  the 
change  in  mean  of  the  residual  is  much  more  significant  than  the  effect  of  the  difference  between  the 
correlation  matrices.  A  similar  argument  can  be  made  for  some  mismodeling  in  the  state  transition 
matrix.  Therefore,  we  are  assuming  that  the  state  estimates  are  large  enough  to  produce  a  significant 
change  in  the  residual  mean  and  that  the  residual  correlation  matrices  are  approximately  equal  when 
we  conduct  a  hypothesis  test  for  a  sensor  failure,  which  is  tantamount  to  assuming  that  the  residual 
mean  effects  will  dominate  the  effects  of  differing  correlation  matrices. 

To  summarize,  we  know  that  the  test  statistic  is  normally  distributed  as  follows; 

*0  :  i(r(t,))  -  N 

D(t.) 

:  L(r(q))  -  N 

We  substitute  this  into  the  known  density  function  for  a  normal  distribution,  which  can  then  be  used 
to  compute  the  threshold  q  for  a  given  Pp^.  Using  the  definition  given  in  Eq  ( 1 10),  we  get: 
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where 


g  = 


r»(t,) 

2 


s 

yig)  =  j 


1 

(2  71 


Y(-.?)=  1  -  Y  (  g)- 


(131) 


Similarly,  we  can  compute  the 


>  =  r _ 1 — 

°  {  (2  7iD(/. 


DW 

.  1 

exp  —  - 

))i/2  2  D(/,.) 


dx 


f 


1- 


0(1/) 


D  (/,)>' 


(2  71 


exp 


=  1-Y 


h  - 


DU/) 


{  O  (/,.)'«  j 


(132) 


For  most  applications,  the  designer  would  set  Pp^  and  Pp  to  achieve  the  desired  hypothesis 
testing  performance,  and  that  would  dictate  a  unique  hypothesis  discrimination  measure,  D  (/, ), 
needed  to  achieve  the  Pp^  and  Pp.  However,  this  does  not  give  us  a  unique  d  (/,)  needed  to  achieve 
the  discrimination  measure,  even  for  a  given  A^'.  Let's  assume  that  one  of  the  hypotheses,  say  h^,  is 
correct,  so  that  /no(^/)  =  T^^is  gives  us  d  (/,)  =  /«,(/,),  which  we  can  use,  along  with  the 
precomputable  to  get  a  class  of  residual  means  that  would  achieve  the  desired  discrimination 
measure.  In  practice,  we  would  choose  a  system  input,  compute  the  residual  mean,  and  compute  the 
discrimination  measure  to  see  it  achieves  the  desired  value.  If  it  doesn't,  we  iterate  on  the  input  until 
the  desired  discrimination  measure  is  achieved. 

The  process  for  computing  the  required  D  (/, )  is  summarized  at  the  end  of  the  next  section. 
Unfortunately,  for  the  flight  control  application  that  we  have  been  considering,  a  single  time  sample 
of  the  residual  will  not  produce  enough  distinguishability  between  two  hypotheses.  We  therefore  will 
look  at  using  multiple  samples  of  the  residual  to  perform  the  hypothesis  testing. 
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simply  to  sum  several  time  samples  of  the  log  likelihood  test  statistic  when  the  hypothesis 
discrimination  is  too  small  to  attain  the  desired  hypothesis  testing  performance.  Therefore,  let 

N-\ 

a(r(r,))  =  Y: 

n-O  (133) 

We  now  need  to  compute  the  mean  of  this  new  test  statistic  under  each  hypothesis.  To  compute  the 
conditional  mean  of  this  test  statistic,  we  will  use  the  conditional  mean  of  each  of  the  terms  in  the 
summation  of  Eq  (133),  conditioned  on  the  measurement  history  up  to  the  last  measurement  before 
the  log  likelihood  test  statistic,  T(r(r..„)),  is  computed.  This  will  allow  us  to  compute  the  log 
likelihood  test  statistic  at  each  time  sample,  and  add  it  to  a  running  sum  of  previous  log  likelihood 
test  statistics,  which  is  precisely  what  we  want  to  implement  in  real  time.  To  denote  the  different 
conditional  means  on  each  of  the  log  likelihood  test  statistics,  we  will  define: 

n  ■  0 

With  this  notation,  we  get: 
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Now  to  compute  the  conditional  covariance  of  this  test  statistic.  As  we  did  above,  we  will 
condition  the  expectation  operator  on  the  measurement  history  up  to  the  latest  data  sample  in  the 
summation  that  we  are  operating  on.  Therefore: 

Since  we  are  developing  a  recursive  relationship  of  9P(r(r,.))  that  is  based  on  the  previous  values,  we 
will  separate  the  sum  as  f^l!''"vs; 


-X  ) 


h!-\ 


n-O 


fxr-i 


1b-0 


=  E 


W-l 


n  -  1 

f 


n  «  1 


(137) 


Expanding  this  relationship  we  get: 

("('/))} 

=  ^Z«,.,)|^'('-('.))'‘  2L(r(/,))E  2^(r(/,.„))» 

fw-i 

f  JV-l 


N- 1 


n  »  1 


n  « 1 


N-  1 


n  =  1 


[^z«i.,){2^(-('.))}] 


rw-i 


f'  (»•(<,-„)) 


^z(^.s){E 


In-  1 
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2  ^z«.,,)^('-(<i))E  ^('■(^i-n))\-Ezct,_,y{Lir(t,))}E,^^  JY:  L(rit,_„)) 


^Z(,.,)  E  -  ^Z(,..)  E  M-(^--n)) 


. .  I 


2 (''('.)>•  E  i(^('.-„)) 

I  n- 1 


If  we  assume  that  the  first  two  terms  of  Eq  (138)  dominate  the  third  term,  so  that: 


In  •  1 


»  2cov2(^  ^JiCrf/,)),  E 


then,  using  this  approximation,  we  get: 


('■('/))}  *  ‘^‘^z(^.r))E  2-(i-(/.  _J) 


D(/.)  »  Lirit._J) 


Since  this  holds  for  any  N,  we  use  this  repeatedly  to  get: 
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cov2(^  ^){S£(r (/.))}*=  D(t.)  *  ,))}.- 

f=  £>(/j)  ♦  25  ( /|.j  )♦•“♦£)(  ,  ) 

N- I 

-  E 

n-O  (141) 

Thus  we  can  approximate  the  covariance  of  the  test  statistic  by  recursively  adding  the  covariances  at 
each  time  sample.  We  will  define  the  disaimination  measure  for  the  test  statistic,  .2(r(t,. )),  as; 


A(h)  =E 

n-0  (142) 

We  have  now  characterized  the  test  statistic,  £([r(i. )),  proposed  in  Eq  (133)  by  defining  the 
mean  of  the  test  statistic,  Eq  (135),  under  both  the  primary  and  alternative  hypotheses,  and  an 
approximation  of  the  covariance  of  the  test  statistic,  Eq  ( 142).  These  equations  were  defined  for  a 
given  number  of  data  samples,  N,  but  what  defines  the  required  number  of  data  samples?  We  will 
now  address  this  question. 

The  performance  of  the  hypothesis  test  is  defined  by  the  probability  of  false  alarms,  Pp^.  and 
the  probability  of  detection,  Pp.  If  we  can  compute  the  inverse  of  the  y  function  defined  in  Eq  ( 1 3 1 ), 
then  we  can  compute  the  test  threshold  and  the  desired  discrimination  measure.  We  will  denote  the 
desired  discrimination  measure  as  Aj,  and  call  it  the  trigger  discrimination  measure  for  reasons 
which  will  become  clear  shortly.  We  will  denote  the  test  threshold  as  tj,  as  before,  and  the  inverse  of 
the  Y  function  as  y'.  We  will  begin  the  computation  of  the  test  threshold  and  trigger  discrimination 
measure  by  using: 


‘  FA 


1-Y 


1/2 


1  -Y(g) 


-  g  =  y  ‘(1  -Pfa) 


(143) 
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and 


Subtracting  Eq  (144)  from  Eq  (143),  we  get: 

^-(^-Ar'^)=  A^'^  =  Y-‘(1-Ppa)-  Y-'(I-Pd) 

"*  ^r  ■  [  Y  (1  -Pfa)  ■  Y  (^  ■  Pd)]  • 
The  test  threshold  can  be  computed  using  the  definition  of  g : 


(144) 


(145) 


(146) 


Some  software  packages  (such  as  MATLAB)  do  not  have  the  inverse  of  the  y  function 
available,  but  they  usually  have  an  inverse  error  function  available.  We  will  extend  the  usual 
development  [53:  1 14-1 15],  Eqs  (143)  through  (146),  to  the  case  where  we  do  not  have  the  inverse 
of  the  Y  function  available  We  can  use  the  following  relations 

\-erf{y)^erfc{y)~~  f  exp  (-t^)dt  «  2  y  (  A  y  ) 

y 

erfmv  (x)  =  y  -  x  =  etf{y)  (147) 

to  compute  the  trigger  discrimination  measure  and  the  test  threshold. 
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First,  we  compute  the  test  threshold  using: 


let  ti 


2y(ii)  =  2y(V^>’)  =  2PpA= 


etf 


1-2? 


FA 


erfinv  ( 1  -  2  Pp^ ) 


JL 


-*  t\  =  'Jl  erfinv  [1-2  Pp^  j . 


(148) 


Then  we  can  compute  the  trigger  discrimination  measure  using: 


let  -*  y  = 


v/5 


-*  2  y(ti  -Ap)  =  2y(v^J')  =  2Pq  =  1 -«;/(>) 


-»  erf 


^-A,^ 


I  V2 


=  1  -  2  Pjj  -»  erfinv  [  1  -  2  Pp  j  = 


n  -  Ap 

~1^ 


-  Ap  =  -n  -v^  er/inv  (1  -2Pp).  (j49) 

Now  we  can  summarize  the  binary  hypothesis  testing  algorithm  that  uses  this  new  test 

statistic. 

Step  1  -  First,  we  would  initialize  the  test  statistic  and  the  discrimination  measure  using: 

a(r(fo))  =  0 

^  ('o)  =  ®  •  (150) 


Step  2  -  Then,  at  each  time  interval  we  would  compute  the  following: 

/«(,(;, )  and  ffi/r, )  defined  in  Eq  (1 19),  using  Eq  (40)  and  Eq  (41), 

</(f, )  and  L(  r(r, ) )  using  Eq  (124)  and  the  uiuiumbered  equation  above  it. 
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D{t,)  =  d{tjA,^d{t,), 

SP(/-(r,))  =  L(r(r,))  +  £P(r(/,.,)). 
and  A(r,.  )  =  D(r,. )  + 


Step  3  -  We  would  then  compare  A(/,. )  to  A^,  as  follows: 

Step  3a  -  If  A(t, )  <  A^,  then  a  hypothesis  test  is  not  performed  for  this  time  iteration. 

Step  3b  -  If  A(t; )  >  A.,-,  then  a  hypothesis  test  is  triggered  (hence  the  name  trigger 
discrimination  measure)  by  comparing  £f(r(t^ ) )  to  p  as  follows: 

Step  3h  -  1  -  If  £fir{i. ) )  >  p,  choose  fi,  (h,  becomes  the  primary  hypothesis)  and  reset 
the  test  statistic  and  discrimination  measure  using: 

<f(r(/.))  =  0 

=  (151) 

Step  3b  -  2  -  If  S£(r{ii ) )  <,  p,  choose  /ig  (/i„  remains  the  primary  hypothesis)  and  reset 
the  test  statistic  and  discrimination  measure  using  Eq  (151). 

This  algorithm  performs  a  binary  hypothesis  test,  but  it  can  be  extended  to  any  A^-ary  test  by 
using  the  algorithm  against  several  alternative  hypotheses.  We  start  by  defining  the  mean  of  the 
residual  given  a  specific  alternative  hypothesis  as: 

(152) 


where  A:  denotes  the  specific  alternative  hypothesis.  We  will  define  S’*  (r(t- ) )  as  the  test  statistic  and 
)  as  the  discrimination  measure  for  the  k""  hypothesis.  Note  that  all  of  the  alternative  hypothesis 
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tests  are  run  in  parallel,  so  all  of  the  test  statistics  and  discrimination  measures  are  computed  in 


parallel  (not  sequentially)  with  those  of  the  other  hypotheses. 

Step  1  -  We  initialize  these  measures  for  each  of  the  alternative  hypotheses: 


=  0 

A,(r,)  =  0 


)  V  Jfc. 


(153) 


Step  2  -  At  each  time  interval  we  compute  /«*(/,),  defined  in  Eq  (152),  using  Eq  (40)  and  Eq  (41),  and 
then  perform  all  of  the  following  computations  for  all  of  the  alternative  hypotheses: 


dt(ti)  =  mj(/.)  -  w,(q) 

=  L^(rCt^))  .  S£*(r(/,.,)) 


(154) 


Initially,  we  will  include  all  of  the  alternative  hypotheses  in  a  list  of  viable  primary 
hypotheses.  We  will  now  start  to  remove  alternative  hypotheses  that  are  not  viable  at  this 
particular  data  sample. 

Step  3  -  Eirst,  we  remove  the  alternative  hypotheses  that  cannot  be  discriminated  from  the  primary 
hypothesis  from  the  list  of  viable  primary  hypotheses.  For  each  of  the  alternative 
hypotheses,  k,  we  would  compare  (t. )  to  A^-,  as  follows: 

Step  3a  -  If  A^  (r, )  <  Ay,  then  a  viable  hypothesis  test  cannot  be  performed  at  this  time 
iteration  for  this  alternative  hypothesis,  so  we  remove  this  alternative  hypothesis  from 
the  list  of  viable  primary  hypotheses. 
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Step  3b  -  If  (t- )  >  Aj,  then  a  viable  hypothesis  test  can  be  performed  at  this  time 
iteration,  so  we  keep  the  alternative  hypothesis  k  on  the  list  of  viable  primary 
hypotheses. 

Step  4  -  Now  we  perform  a  hypothesis  test  between  the  current  primary  hypothesis  and  all  of  the 

alternative  hypotheses  left  on  the  viable  primary  hypothesis  list.  For  each  of  the  alternative 
hypotheses  left  on  the  list  of  viable  primary  hypotheses,  we  perform  a  hypothesis  test  by 
comparing  (r(t. ) )  to  rj  as  follows: 

Step  4a  -  If  a’t  (r(/, ) )  <  rj ,  the  alternative  hypothesis  is  removed  from  the  list  of  viable 
primary  hypotheses  and  its  test  statistic  and  discrimination  measure  are  reset,  only  for 
this  particular  alternative  hypothesis.  Note  that  the  alternative  hypotheses  that  were 
removed  from  the  list  of  viable  primary  hypotheses  did  not  trigger  a  hypothesis  test,  so 
we  do  noi  reset  their  test  statistic  and  discrimination  measure. 

Step  4b  -  If  (r(/, ) )  >  r| ,  keep  the  alternative  hypothesis  on  the  list  of  viable  primary 
hypotheses. 

Step  5  -  We  now  have  two  possible  cases,  either  the  list  of  viable  primary  hypotheses  still  has  one  or 
more  entries,  or  it  does  not. 

Step  5a  -  If  the  list  of  viable  primary  hypotheses  is  empty,  then  the  current  primary 
hypothesis  is  still  in  effect. 

■Step.  5b  -  If  the  list  of  viable  primary  hypotheses  is  not  empty,  then  we  choose  the  hypothesis 
that  has  the  largest  test  statistic  and  reset  all  of  the  test  statistics  and  discrimination 
measures  using; 


0 

0 


V  k. 


(155) 
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Step  4a  allows  the  other  alternative  hypotheses  to  trigger  the  hypothesis  test  at  the 
appropriate  number  of  iterations  to  provide  the  desired  hypothesis  testing  performance.  Thus,  the 
hypotheses  that  have  discrimination  measures  that  grow  faster  because  they  are  easily  distinguished 
from  the  primary  hypothesis  would  be  tested  before  the  less  easily  distinguished  hypotheses. 

For  example,  let  the  primary  hypothesis  be  a  no  flight  control  failure  and  two  alternative 
hypotheses  a  left  and  right  aileron  failure.  If  the  right  aileron  failure  hypothesis  discrimination 
measure  grew  faster  than  the  left  aileron  failure  hypothesis  disaimination,  a  no  failure  hypothesis 
versus  a  right  aileron  failure  hypothesis  test  would  be  triggered  before  the  no  failure  vs.  left  aileron 
failure  hypothesis  test.  Once  the  test  is  triggered,  if  the  no  failure  hypothesis  is  chosen,  then  the  test 
discrimination  measure  for  the  right  aileron  failure  hypothesis  is  reset.  The  test  discrimination 
measure  for  the  left  aileron  failure  hypothesis  is  not  reset  so  it  will  continue  to  grow  until  a  no  failure 
vs.  left  aileron  failure  hypothesis  test  is  triggered.  Note  that  if  the  left  aileron  failure  hypothesis 
discrimination  measure  were  reset  in  this  case,  a  no  failure  vs.  left  aileron  failure  hypothesis  test 
would  never  be  triggered.  If  the  no  failure  vs.  right  aileron  failure  hypothesis  test  is  triggered  and  the 
alternative  hypothesis  (right  aileron  failure)  is  chosen,  then  aj]  of  the  test  measures  and  hypothesis 
discrimination  measures  are  reset  because  we  would  have  a  new  primary  hypothesis  against  which  we 
will  test  the  alternative  hypotheses. 

Note  that  this  is  not  truly,  an  optimal  failure  hypothesis  testing  structure  because  we  are  only 
testing  the  primary  hypothesis  against  the  alternative  hypotheses,  and  we  are  not  testing  any 
alternative  hypotheses  against  each  other.  In  the  example  above,  we  are  not  testing  the  left  vs.  right 
aileron  failures.  This  could  artificially  boost  the  likelihood  of  an  incorrect  hypothesis  in  the  event  of 
a  missed  detection.  For  example,  if  a  no-failure  vs.  right  aileron  hypothesis  test  is  triggered,  and  a 
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right  aileron  failure  is  mistakenly  missed,  then  the  left  aileron  failure  discrimination  measure  is  not 
reset  and  most  likely  will  reach  its  discrimination  trigger  before  the  right  aileron  failure 
discrimination  measure  (which  was  reset  after  the  test  was  conducted)  reaches  its  trigger.  If  there 
happens  to  be  ambiguity  between  the  two  hypotheses  (right  vs.  left  aileron  failure),  the  right  aileron 
vs.  no-failure  hypothesis  test  will  likely  choose  the  incorrect  hypothesis  (the  right  aileron  failure). 
This  occurs  because  we  are  not  tracking  the  discrimination  between  the  alternative  hypotheses,  which 
would  require  )*((/f-l)  (where  K  is  the  number  of  hypotheses)  different  hypothesis  triggers, 
discrimination  measures,  and  tests.  If  it  is  known  that  certain  alternate  hypotheses  have  ambiguity 
problems,  the  probability  of  missed  detection  (l-P^)  for  those  hypotheses  could  be  set  lower  than  the 
other  hypotheses  to  decrease  the  probability  of  this  type  of  problem. 

The  timing  of  the  failure  could  influence  the  amount  of  time  needed  to  identify  the  failure. 
Ideally,  if  the  failure  occurs  right  when  the  test  disaiminalion  measure  is  started,  this  algorithm  will 
trigger  a  hypothesis  test  just  at  the  time  when  the  hypothesis  is  distinguishable  from  the  primary 
hypothesis,  say  M  time  samples  after  the  failure.  However,  if  the  failure  occurred  after  the  lest 
discrimination  measure  has  already  started,  the  discrimination  measure  will  mgger  a  test  before  the 
two  hypotheses  are  truly  distinguishable,  so  it  is  quite  possible  for  the  algorithm  to  miss  the  failure 
detection  at  this  first  test,  but  it  should  pick  up  the  failure  M  lime  samples  later.  Thus,  the  worst  case 
delay  in  the  failure  detection  of  this  algorithm  would  be  2M-\. 

3.4  Chapter  Summary 

We  have  established  the  theoretical  and  algorithmic  contributions  of  this  research  in  this 
chapter.  We  started  by  presenting  the  Kalman  filter  equations  and  laid  out  the  nomenclature  that  we 
use  to  represent  the  mismodeling  in  the  models  that  are  used  to  develop  the  Kalman  filter  equations. 
The  specific  mismodeling  cases  that  we  represented  are  a  mismodeled  input  matrix,  output  matrix. 
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and  state  transition  matrix.  We  developed  expressions  for  the  mean  and  covariance  of  the  residual 
for  any  of  these  three  types  of  mismodeling,  and  showed  how  the  expressions  reduced  for  specific 
cases  of  mismodeling.  In  Section  3.2.5  we  developed  an  equivalent  structure  to  the  Standard  Kalman 
Filter  Bank  (SKFB)  that  uses  the  residual  and  state  estimates  from  a  single  Kalman  filter,  to  compute 
the  equivalent  residuals  of  all  other  Kalman  filters  with  known  model  differences,  which  would 
usually  be  implemented  in  an  MMAE  algorithm.  We  developed  a  Residual  Correlation  Kalman 
Filter  Bank  (RCKFB)  in  Section  3.2.6,  that  estimates  the  power  spectral  density  of  the  residuals  and 
feeds  this  to  the  hypothesis  testing  algorithm,  motivated  by  the  idea  that  this  structure  would  exploit 
the  known  characteristics  of  the  correlation  of  the  residual  to  identify  failures.  Section  3.3. 1  defined 
the  Standard  Hypothesis  Testing  Algorithm  (SHTA)  and  described  several  of  the  modifications  that 
have  been  researched.  We  developed  a  Neyman-Pearson  based  Hypothesis  Testing  Algorithm 
(NPHTA)  in  Section  3.3.2,  where  we  extended  the  usual  single  time  sample,  binary  hypothesis  test 
formulation,  to  multiple  time  sample,  multiple  hypothesis  testing.  These  developments  suggest  an 
algorithm  architecture  of  an  "MM  AF-like"  algorithm,  but  with  the  many  elemental  filters  and 
residuals  of  a  usual  MMAE  replaced  by  one  filler,  linear  transformations  (optimized  "matched 
filters")  and  equivalent  residuals.  An  eventually  implemented  system  might  be  composed  of  an 
algorithm  structure  based  on  wither  the  SKFB  or  RCKFB  (or  both),  possibly  using  equivalent 
residuals  rather  than  actual  residuals,  and  perhaps  using  the  NPHTA  rather  than  the  SHTA.  The 
next  chapter  applies  these  various  algoritlims  to  a  specific  sensor/actuator  failure  detection  problem 
to  assess  relative  strengths  of  the  various  configurations. 
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IV.  Results 


4. 1  Chapter  Overview 

In  this  chapter  we  present  the  results  of  our  research.  In  Section  4.2,  we  verify  our 
development  of  an  expression  for  the  mean  of  the  residual  (Eq  (40)  and  Eq  (41)  from  Section  3.2.4) 
and  the  covariance  of  the  residual  (Eq  (34)  from  Section  3.2.3).  We  present  a  graphical  comparison 
of  the  time  history  of  the  predicted  mean  and  covariance  for  each  element  of  a  residual  along  with  a 
single  sample  of  a  residual  from  a  Kalman  filter.  A  more  detailed  comparison  of  the  mean  and 
covariance  of  the  residual  is  made  by  temporally  averaging  five  samples  of  the  residual  and  then 
ensemble  averaging  these  temporal  averages,  and  ensemble  averaging  10  samples  of  the  residual  and 
temporally  averaging  these  ensemble  (sample)  averages.  These  empirically  computed  statistics  are 
compared  to  the  predicted  values  of  the  mean  and  covariance.  This  comparison  is  made  for  the 
specific  cases  of  no  mismodeling  in  the  Kalman  filler  model,  a  mismodeled  input  control  matrix 
(actuator  failure),  and  a  mismodeled  output  matrix  (sensor  failure).  Time  limitations  prevented  the 
computation  of  results  for  the  mismodeled  state  transition  matrix  test  case.  The  development  of 
equivalent  residuals,  presented  in  Section  3.2.5,  is  verified  in  Section  4.3.  We  compare  the  residuals 
from  a  bank  of  several  Kalman  filters  with  the  equivalent  residuals  that  are  computed  using  a  single 
Kalman  filter  residual  and  known  model  differences  between  the  Kalman  filters.  The  comparison  is 
made  for  the  same  specific  cases  of  mismodeling  that  were  presented  in  Section  4.2.  These  results 
suggest  a  new  filter  bank  architecture  that  replaces  the  many  filters  and  residuals  with  a  single 
Kalman  filter,  several  linear  transforms,  and  equivalent  residuals. 
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Section  4.4  is  a  comparison  of  the  failure  identification  performance  of  the  three  different 
MMAE  structures  that  were  developed:  Standard  Kalman  Filter  Bank  (SKFB)  with  a  Standard 
Hypothesis  Testing  Algorithm  (SHTA),  Residual  Correlation  Kalman  Filter  Bank  (RCKFB)  with  a 
SHTA,  and  SKFB  with  a  Neyman-Pearson  Hypothesis  Testing  Algorithm  (NPHTA).  Time 
limitations  prevented  the  development  and  testing  of  a  fourth  configuration,  namely  a  RCKFB  with 
NPHTA.  The  testing  of  the  failure  identification  performance  was  accomplished  for  various  actuator 
failures,  which  were  modeled  hy  a  cnhimn  of  (he  conirol  input  matrix  being  set  to  zero.  Previous 
research  [19,  31, 41]  has  shown  that  these  were  the  most  difficult  failure  modes  to  identify  for  this 
particular  aircraft,  so  we  chose  to  test  the  performance  for  these  failures  to  be  able  to  characterize  the 
worst  case  performance  of  these  various  structures.  The  comparison  of  these  different  architectures 
suggests  an  algorithm  structure  that  uses  both  the  SKFB  and  RCKFB,  perhaps  with  equivalent 
residuals  instead  of  the  actual  residuals,  using  a  NPHTA  for  the  best  failure  identification 
performance. 

4.2  Characterization  of  the  Re.sidual 

In  Section  3.3,  we  developed  a  recursive  algorithm  for  estimating  the  mean  and  covariance  of 
the  Kalman  filter  residuals  under  various  forms  of  hypothesized  mismodeling.  We  present  a 
comparison  of  these  estimates  with  actual  Kalman  filter  residuals.  The  comparisons  are  presented 
for  three  hypotheses;  no  mismodeling,  an  input  control  matrix  mismodeling  {AB  *  0),  and  an  output 
matrix  mismodeling  (AH  *  0).  For  the  flight  control  failure  identification  application,  these 
hypotheses  correspond  to  no  flight  control  failure,  an  actuator  failure,  and  a  sensor  failure, 
respectively. 

The  application  that  we  used  to  test  the  theory  developed  in  this  dissertation  consisted  of  a 
bank  of  16  Kalman  filters,  each  modeling  a  different  hypothesis  (1  no  failure,  6  actuator  failures,  and 
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9  sensor  failures).  The  complete  set  of  results  for  this  section  would  require  us  to  present  256  test 
cases,  16  Kalman  filter  residuals  for  each  of  the  16  hypotheses.  Therefore,  we  will  present 
representative  results  to  show  the  accuracy  of  the  estimates  of  a  Kalman  filter  residual  for  the  three 
cases  of  mismodeling. 

The  plots  presented  below  show  each  of  the  scalar  elements  of  a  single  sample  of  the  9- 
element  residual  in  a  separate  plot.  The  elements  are  shown  as  a  function  of  time,  with  the 
simulation  time  being  8  seconds.  The  covariance  of  the  residual  is  precomputed  usir.g  .Section  .T2.3. 
For  each  element  of  the  residual,  wc  superimpose  on  the  plot  of  the  single  sample  of  the  residual 
process,  the  computed  mean,  the  computed  mean  plus  the  square  root  of  the  computed  variance 
(corresponding  to  one  computed  standard  deviation  for  the  residual  element),  and  the  computed  mean 
minus  one  standard  deviation.  Since  the  residual  is  assumed  to  be  Gaussian,  we  would  expect  that 
about  68%  of  the  time  samples  of  the  residual  would  lie  within  the  plotted  mean  ±  one  standard 
deviation  bounds.  Thus,  we  can  visually  compare  the  residual  with  its  computed  mean  and 
covariance. 

A  more  detailed  analysis  of  the  residual  characteristics  is  accomplished  by  comparing 
various  empirical  statistics  to  the  computed  mean  and  covariance.  First,  we  computed  the  temporal 
average  (ergodic  approximation  of  the  mean)  and  the  temporally-computed  standard  deviation  of 
each  of  five  realizations  of  the  residual  process.  Then  we  averaged  these  temporal  averages  to  get  the 
ensemble-temporal  average  of  the  mean  and  standard  deviation.  Next,  we  computed  the  ensemble 
(sample)  average  of  10  realizations  of  the  residual  process,  which  gave  us  a  time  history  average  of 
both  the  mean  and  standard  deviation.  These  ensemble  (sample)  averages  were  averaged  temporally 
to  obtain  the  temporal-ensemble  averages.  These  two  types  of  averages  are  compared  to  the 
predicted  mean  and  standard  deviation. 
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We  chose  to  examine  the  residual  vector  from  the  Kalman  filter  that  models  a  fully 
functional  aircraft  (no  flight  control  failures).  The  true  system  model  assumes  a  fully  functional 
airaaft  for  the  first  second  of  simulation,  followed  by  a  specific  single  flight  control  failure.  Thus, 
for  the  first  second  of  simulation,  the  chosen  Kalman  filter  has  no  mismodeling,  after  which  a 
specific  mismodeling  occurs,  which  will  be  the  difference  between  the  fully  functional  aircraft  model 
and  the  actual  flight  control  failure  model.  We  will  look  at  specific  cases  of  mismodeling  in  the 
following  sections. 

4.2.1  No  Mismodeling  (No  Failure).  We  start  by  comparing  the  residual  from 

the  Kalman  filter  that  models  a  fully  functional  aircraft  (no  flight  control  actuator  or  sensor  failures) 
where  the  model  is  accurate  throughout  the  simulation  time.  The  mean  and  covariance  were 
computed  using  Eq  (40)  and  Eq  (41)  with  AB=:0,  A//=0,  and  A$=  0.  The  results  for  this  specific 
case  of  mismodeling  were  developed  in  Section  3.2.4. 1,  where  we  found  that  the  residual  should  be 
zero-mean  for  all  time.  Figure  12  shows  that  the  residual  is  zero-mean  throughout  the  simulation  and 
approximately  68%  of  the  time  samples  are  within  the  computed  standard  deviation  bounds. 

To  facilitate  a  more  accurate  comparison  of  the  residual  with  its  computed  characteristics, 
we  present,  in  Table  1,  the  computed  the  mean  and  standard  deviation  of  the  scalar  components,  the 
ensemble-temporal  averages,  and  the  temporal-ensemble  averages.  The  computed  standard  deviation 
of  each  element  was  calculated  by  computing  the  diagonal  elements  of  the  matrix  square  root  of  the 
computed  covariance  matrix  (which  was  shown  to  be  time  invariant). 

These  results  clearly  indicate  that  the  residual  was  well  characterized  for  the  case  of  no 
mismodeling  in  the  Kalman  filter  model.  Note  that  these  averages  (computed  using  both  methods) 
may  differ  because  of  the  different  number  samples  used  for  each.  Also,  these  averages  are 
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Fully  Functional  Kalman  Filter  with  No  Mismodeling 
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Table  1 


Kalman  Filter  Residual  Statistical  Comparison  when  No  Mismodeling  Exists. 


Residual 

Element 

Computed 

Average 

Ensemble 

Temporal 

Average 

Temporal 

Ensemble 

Average 

Computed 

0 

Ensemble 

Temporal 

0 

Temporal 

Ensemble 

0 

Forward  Velocity 

0 

-1.34  E-1 

-2.24  E-1 

3.00 

2.97 

2.96 

Angle  of  Attack 

0 

3.70  E-3 

-8.04  E-4 

2.96  E-2 

2.95  E-2 

2.95  E-2 

Pitch  Rate 

0 

4.09  E-4 

4.65  E-3 

3.55  E-2 

3.53  E-2 

3.50  E-2 

Pitch  Angle 

0 

5.76  E-4 

3.80  E-5 

2.92  E-3 

2.87  E-3 

2.97  E-3 

Sideslip  Angle 

0 

-3.64  E-3 

-6.51  E-3 

1.44  E-1 

1.44  E-1 

1.42  E-1 

Roll  Rate 

0 

-2.30  E-4 

-1.75  E-3 

3.46  E-2 

3.41  E-2 

3.48  E-2 

Roll  Angle 

0 

-1.62  E-5 

-7.12  E-5 

3.94  E-3 

3.80  E-3 

3.99  E-3 

Yaw  Rate 

0 

9.67  E-5 

5.66  E-4 

5.32  E-2 

5.57  E-2 

5.40  E-2 

Yaw  Angle 

0 

-5.81  E-5 

-5.34  E-5 

3.22  E-3 

3.12  E-3 

3.20  E-3 

dominated  significantly  by  the  standard  deviation  values,  as  we  would  expect  from  the  computed 
mean  (i.e.,  zero)  and  standard  deviations. 

4.2.2  Mismodeled  Input  Control  Matrix  (Actuator  Failure).  Wenow 

compare  the  residual  from  the  Kalman  filter  that  uses  a  fully  functional  aircraft  model  to  the 
computed  mean  and  covariance  assuming  that  the  control  input  matrix  is  mismodeled  starting  at  1 
second  into  the  simulation.  This  situation  corresponds  to  an  actuator  failure  occurring  at  1  second. 
We  will  look  at  an  elevator  actuator  failure,  an  aileron  actuator  failure,  and  a  rudder  actuator  failure. 
Similar  results  were  obtained  using  the  Kalman  filter  that  incorrectly  modeled  the  failure  status  from 
the  start  of  simulation  to  1  second,  and  then  correctly  modeled  the  failure  status  from  then  to  the  end 
of  the  simulation.  An  example  of  this  case,  where  the  filler  model  is  based  on  a  left  elevator  failure 
and  the  left  elevator  failure  occurs  at  1  second,  is  shown  later  in  Figure  21,  which  shows  that  the 


107 


residual  initially  is  not  zero-mean,  but  eventually,  after  the  failure  occurs,  becomes  zero-mean.  For 
all  of  these  cases,  the  mean  of  the  residual  was  computed  using  Eq  (40)  and  Eq  (41)  and  the  expected 
covariance  was  developed  in  Section  3.2.3.  In  Section  3. 2.4.2  we  found  that  AB^O  results  in  a 
residual  mean  that  has  components  of  the  control  input,  including  purposeful  dithers  that  are  put  into 
the  input  to  enhance  failure  identification.  The  inputs  that  were  used  for  this  research  are  shown  in 
Figure  13. 

The  results  for  the  elevator  failure  case,  are  graphically  shown  in  Figure  14.  For  this  case  we 
have  no  mismodeling  up  to  1  second  in  the  simulation,  then  mismodeling  occurs  in  the  column  of  the 
control  input  matrix  that  corresponds  to  an  elevator,  a  left  elevator  in  this  case.  The  left  elevator 
input  is  quite  evident  in  the  residual,  particularly  in  the  pitch  rate  and  pitch  angle  where  there  is  a 
significant  change  in  scale.  Also  the  forward  velocity  and  angle  of  attack  show  a  slight  change  in 
mean  corresponding  to  the  left  elevator  input.  The  other  elements  of  the  residual  do  not  show  the  left 
elevator  input  because  the  aircraft  model  assumes  there  is  no  cross-coupling  between  the  pitch  axis 
and  the  other  (lateral-directional)  axes. 

'  These  are  precisely  the  results  that  we  expected.  If  an  elevator  actuator  failed,  then  the 
aircraft  will  not  respond  to  the  elevator  input.  The  Kalman  filter  with  the  model  based  on  a  fully 
functional  aircraft,  would  predict  the  measurements  based  on  the  aircraft  responding  to  the  elevator 
input.  Since  the  aircraft  is  nM  responding  to  this  input,  the  difference  between  the  Kalman  filter's 
prediction  of  the  measurements  and  the  actual  measurements  (the  residual)  would  show  the  aircraft 
response.  For  this  particular  case,  a  sinusoidal  elevator  input  would  case  a  change  in  the  pitch  rate, 
pitch  angle,  angle  of  attack,  and  forward  velocity.  Since  the  fully  functional  Kalman  filter  would 
predict  these  responses,  but  the  aircraft  is  actually  not  responding  to  this  elevator  input,  the 
difference  between  the  filter's  prediction  and  the  actual  measurements  would  show  this  expected 
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Figure  13.  Control  Inputs,  Plotted  in  Radians. 
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Figure  14.  Residual  and  Computed  Residual  Mean  and  Standard  Deviation  of  the 
Fully  Functional  Kalman  Filter  with  a  Mismodeled  Elevator  Input. 


no 


response  to  the  elevator  input.  This  is  precisely  what  we  found  in  Figure  14.  Similar  results  were 
obtained  for  a  right  elevator  failure. 

To  facilitate  further  analysis,  we  computed  the  bias-compensated  residual  by  differencing  the 
computed  mean  and  the  value  of  the  realization  of  the  residual  process  for  each  time  sample: 

?(',)  =  »«c(h)  - '■(h)  (156) 


where  is  the  computed  residual  mean.  If  the  computed  residual  mean  accurately  represents  the 
residual  mean,  then  we  would  expect  that  the  ensemble  mean  of  the  temporal  average  of  this 
difference  to  be: 
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Likewise,  the  temporal  average  of  the  ensemble  mean  of  this  difference  is: 
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Since  the  only  random  part  of  5  is  the  residual,  the  covariance  of  this  bias-compensated  residual  is 
the  covariance  of  the  residual. 


Ill 


Thus,  if  we  compute  the  temporal-ensemble  average  (compute  the  ensemble  average  first 
and  then  compute  the  temporal  average  of  the  ensemble  average)  and  the  ensemble-temporal  average 
(compute  the  temporal  average  of  each  realization  and  then  compute  the  ensemble  average  of  these 
temporal  averages)  we  would  expect  the  same  results  that  were  predicted  for  No  Mismodeling  case, 
shown  in  Table  1. 

We  computed  the  results  for  three  cases,  an  elevator  failure,  an  aileron  failure,  and  a  rudder 
failure.  The  results  for  the  first  case,  an  elevator  failure,  are  shown  in  Table  2.  Clearly,  the  residual 
is  accurately  characterized. 

Table  2 


Kalman  Filter  Residual  Statistical  Comparison  for  the  Case  of  an  Elevator  Failure. 


Residual 

Element 

Computed 
Average 
of  ^ 

Ensemble 

Temporal 

Average 

of^ 

Temporal 

Ensemble 

Average 

of^ 

Computed 

0  of  ^ 

Ensemble 

Temporal 

0  of  ^ 

Temporal 

Ensemble 

0  of  ^ 

Forward  Velocity 

0 

-2.05  E-1 

-2.35  E-1 

3.00 

3.00 

2.98 

Angle  of  Attack 

0 

-4.83  E-3 

-1.99  E-4 

2.96  E-2 

2.92  E-2 

2.96  E-2 

Pitch  Rate 

0 

-3,70  E-3 

-1.30  E-3 

3.55  E-2 

3.54  E-2 

3.60  E-2 

Pitch  Angle 

0 

5.76  E-4 

-9.60  E-5 

2.92  E-3 

2.76  E-3 

2.87  E-3 

Sideslip  Angle 

0 

3.94  E-3 

2.45  E-4 

1.44  E-1 

1.43  E-1 

1.42  E-1 

Roll  Rate 

0 

-2.30  E-4 

-4.03  E-4 

3.46  E-2 

3.38  E-2 

3.44  E-2 

Roll  Angle 

0 

3.55  E-4 

-7.78  E-5 

3.94  E-3 

3.87  E-3 

3.97  E-3 

Yaw  Rate 

0 

3.22  E-5 

5.66  E-4 

5.32  E-2 

5.49  E-2 

5.43  E-2 

Yaw  Angle 

0 

-5.81  E-5 

8.08  E-6 

3.22  E-3 

3.15  E-3 

3.18  E-3 
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The  next  case  corresponds  to  an  aileron  failure.  Again,  we  have  no  mismodeling  up  to  1 
second  in  the  simulation,  then  mismodeling  occurs  in  the  column  of  the  control  input  matrix  that 
corresponds  to  an  aileron.  In  Figure  15  we  graphically  show  the  results;  the  input  to  the  left  aileron 
is  clearly  evident  in  the  residual,  particularly  in  the  roll  rate  and  roll  angle,  and  to  a  lesser  degree  in 
the  yaw  rate  and  yaw  angle.  Note  the  dramatic  change  in  scale  for  these  elements  of  the  residual 
vector.  Just  like  in  the  elevator  failure  case,  the  aircraft  is  not  responding  to  the  left  aileron  input,  but 
the  fully  functional  Kalman  filter  is  predicting  that  the  aircraft  will  respond  to  this  input.  Thus,  the 
difference  between  the  actual  measurements  (which  show  the  aircraft's  lack  of  response)  and  the 
predicted  measurements  (which  assume  that  the  aircraft  is  responding)  will  reflect  the  aircraft 
response  to  a  left  aileron  input.  The  expected  aircraft  response  to  a  sinusoidal  left  aileron  input 
would  be  a  sinusoidal  change  in  the  roll  rate  and  roll  angle.  Also,  since  this  is  an  input  to  only  one 
aileron,  we  expect  a  small  sinusoidal  change  in  the  yaw  rate  and  yaw  angle.  This  is  precisely  what 
was  observed  in  Figure  15.  Similar  results  were  obtained  for  a  right  aileron  failure. 

We  also  computed  and  tabulated  the  ensemble-temporal  average  and  temporal-ensemble 
average  of  the  bias-compensated  residual.  These  results,  shown  in  Table  3,  show  that  the 
characteristics  of  the  residual  are  accurately  computed  for  the  case  of  an  aileron  failure. 

The  last  case  corresponds  to  a  rudder  failure.  Again,  we  have  no  mismodeling  up  to  1 
second  in  the  simulation,  then  mismodeling  occurs  in  the  column  of  the  control  input  matrix  that 
corresponds  to  a  rudder.  In  Figure  16  we  show  the  results  for  a  left  rudder  failure.  The  input  to  the 
left  rudder  is  clearly  evident  in  the  residual,  particularly  in  the  yaw  rate  and  yaw  angle  with  some 
cross-coupling  to  the  roll  rate  and  roll  angle.  Again,  note  the  change  in  scale  for  these 
residualelements.  This  aircraft  has  two  rudders  (similar  to  the  A- 10  aircraft),  so  if  one  rudder  is  not 
responding  to  a  sinusoidal  input,  the  other  rudder  would  cause  not  only  sinusoidal  changes  in  yaw 
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Figure  15.  Residual  and  Computed  Residual  Mean  and  Standard  Deviation  of  the 
Fully  Functional  Kalman  Filter  with  a  Mismodeled  Aileron  Input. 
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Tables 


Kalman  Filter  Residual  Statistical  Comparison  for  the  Case  of  an  Aileron  Failure. 


Residual 

Element 

Computed 

Average 

of^ 

Ensemble 
Temporal 
Average 
of  5 

Temporal 

Ensemble 

Average 

ofC 

Computed 

0  of  5 

Ensemble 

Temporal 

0  of  5 

Temporal 
Ensemble 
o  of  ^ 

Forward  Velocity 

0 

2.24  E-1 

-9.50  E-2 

3.00 

2.96 

2.96 

Angle  of  Attack 

0 

-1.76  E-3 

-1.08  E-4 

2.96  E-2 

2.91  E-2 

2.78  E-2 

Pitch  Rate 

0 

1.38  E-3 

6.54  E-4 

3.55  E-2 

3.58  E-2 

3.44  E-2 

Pitch  Angle 

0 

4.68  E-5 

3.23  E-4 

2.92  E-3 

2.85  E-3 

2.89  E-3 

Sideslip  Angle 

0 

1.59  E-2 

1.46  E-3 

1.44  E-1 

1.42  E-1 

1.42  E-1 

Roll  Rate 

0 

1.10  E-3 

9.89  E-4 

3.46  E-2 

3.45  E-2 

3.48  E-2 

Roll  Angle 

0 

-2.76  E-4 

2.05  E-3 

3.94  E-3 

3.67  E-3 

3.75  E-3 

Yaw  Rate 

0 

1.13  E-3 

2.86  E-3 

5.32  E-2 

5.52  E-2 

5.46  E-2 

Yaw  Angle 

0 

3.18  E-3 

6.09  E-4 

3.22  E-3 

3.19  E-3 

3.19  E-3 

Table  4 

Kalman  Filter  Residual  Statistical  Comparison  for  the  Case  of  a  Rudder  Failure. 


Residual 

Element 

Computed 

Average 

of? 

Ensemble 

Temporal 

Average 

of? 

Temporal 

Ensemble 

Average 

of?'' 

Computed 
o  of  ? 

Ensemble 
Temporal 
o  of  ? 

Temporal 
Ensemble 
o  of  ? 

Forward  Velocity 

0 

-2.01  E-1 

-1.85  E-1 

3.00 

2.95 

2.98 

Angle  of  Attack 

0 

-1.12  E-3 

5.60  E-4 

2.96  E-2 

2.97  E-2 

2.98  E-2 

Pitch  Rate 

0 

2.61  E-3 

2.66  E-3 

3.55  E-2 

3.59  E-2 

3.53  E-2 

Pitch  Angle 

0 

-7.83  E-5 

4.19  E-5 

2.92  E-3 

2.94  E-3 

2.86  E-3 

Sideslip  Angle 

0 

-2.16  E-3 

-3.79  E-4 

1.44  E-1 

1.43  E-1 

1.43  E-1 

Roll  Rate 

0 

6.07  E-4 

6.09  E-4 

3.46  E-2 

3.53  E-2 

3.48  E-2 

Roll  Angle 

0 

1.20  E-3 

3.16  E-4 

3.94  E-3 

3.51  E-3 

3.89  E-3 

Yaw  Rate 

0 

3.01  E-3 

2.09  E-3 

5.32  E-2 

5.35  E-2 

5.28  E-2 

Yaw  Angle 

0 

5.71  E-3 

1.95  E-4 

3.22  E-3 

3.15  E-3 

3.27  E-3 
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Figure  16.  Residual  and  Computed  Residual  Mean  and  Standard  Deviation  of  the 
Fully  Functional  Kalman  Filter  with  a  Mismodeled  Rudder  Input. 
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rate  and  yaw  angle,  but  some  small  cross-coupling  would  occur  in  roll  rate  and  roll  angle.  Again, 
since  the  fully  functional  Kalman  filter  predicts  that  the  airaaft  is  responding  to  the  input,  but  the 
airaaft  actually  is  not,  the  residual  would  reflect  the  airaaft  response  to  a  single  rudder  input.  This 
is  precisely  what  is  reflected  in  Figure  16.  As  before,  the  difference  results  were  computed  and  are 
tabulated  in  Table  4.  Clearly,  the  residual  is  accurately  characterized.  Similar  results  were  obtained 
for  a  right  rudder  failure. 

4,2.3  Mismodeled  Output  Matrix  (Sen.sor  Failure).  We  now  compare  the 

residual  from  the  Kalman  filter  that  uses  a  fully  functional  aircraft  model  to  the  computed  mean  and 
covariance,  assuming  that  the  output  matrix  is  mismodeled  starting  at  1  second  into  the  simulation. 
This  situation  corresponds  to  a  sensor  failure  occurring  at  1  second.  We  look  separately  at  a  failure 
in  the  forward  velocity  sensor  and  the  pitch  rate  sensor.  Once  again,  the  mean  of  the  residual  was 
computed  using  Eq  (40)  and  Eq  (41)  and  the  expected  covariance  was  developed  in  Section  3.2.3. 
The  failure  mode  that  is  represented  by  this  mismodeling  in  the  output  matrix  (a  row  of  the  matrix  is 
zero)  corresponds  to  a  "fail  to  nominal,"  rather  than  "fail  to  zero,"  since  the  Kalman  filter  model  is  a 
perturbation  model.  The  mismodeling  in  the  output  matrix  zeros  out  any  non-zero  value  of  the 
perturbation  state,  so  any  perturbations  from  the  nominal  are  set  to  zero  for  this  failure  mode.  The 
results  for  a  forward  velocity  sensor  failure  are  graphically  shown  in  Figure  17,  where  the  forward 
velocity  residual  element  is  clearly  no  longer  zero-mean  (note  the  dramatic  change  in  scale).  The 
aircraft  model  did  not  include  a  thrust  input,  so  we  consistently  noticed  that  the  forward  velocity 
decreased  throughout  the  simulation  (since  there  was  no  thrust  to  counteract  the  drag  caused  by  the 
actuator  inputs).  Once  the  forward  velocity  sensor  failed,  this  decrease  in  velocity  was  not  measured, 
but  the  fully  functional  Kalman  filter  continued  to  predict  that  the  forward  velocity  would  decrease. 
Therefore,  the  difference  between  the  actual  measurement  and  the  filter's  prediction  of  that 
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Fully  Functional  Kalman  Filter  with  a  Mismodeled  Forward  Velocity  Sensor. 
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Table  5 


Kalman  Filter  Residual  Statistical  Comparison  for  the  Case  of  a  Forward  Velocity  Sensor  Failure. 


Residual 

Element 

Computed 

Average 

of? 

Ensemble 

Temporal 

Average 

of? 

Temporal 

Ensemble 

Average 

of? 

Computed 
o  of  ? 

Ensemble 
Temporal 
o  of  ? 

Temporal 
Ensemble 
o  of  ? 

Forward  Velocity 

0 

-2.01  E-1 

2.27  E-1 

3.003 

2.98 

3.05 

Angle  of  Attack 

0 

-3.79  E-4 

-1.01  E-3 

2.96  E-2 

2.89  E-2 

2.97  E-2 

Pitch  Rate 

0 

8.55  E-4 

9.41  E-2 

3.55  E-2 

3.47  E-2 

3.66  E-2 

Pitch  Angle 

0 

-1.01  E-5 

-6.13  E-3 

2.92  E-3 

2.87  E-3 

2.94  E-3 

Sideslip  Angle 

0 

6.82  E-4 

1.64  E-3 

1.44  E-1 

1.43  E-1 

1.41  E-1 

Roll  Rate 

0 

1.55  E-4 

4.86  E-4 

3.46  E-2 

3.41  E-2 

3.38  E-2 

Roll  Angle 

0 

2.62  E-4 

1.76  E-4 

3.94  E-3 

3.86  E-3 

3.97  E-3 

Yaw  Rate 

0 

-9.14  E-4 

-1.73  E-3 

5.32  E-2 

5.36  E-2 

5.29  E-2 

Yaw  Angle 

0 

2.87  E-5 

1.12  E-4 

3.22  E-3 

3.20  E-3 

3.18  E-3 

measurement  is  clearly  reflected  in  the  residual.  The  difference  between  the  computed  mean  and  the 
actual  residual  was  computed  and  the  temporal  mean  and  standard  deviation  were  calculated.  The 
results  are  tabulated  in  Table  5.  Once  again,  the  characteristics  of  die  residual  are  accurately 
computed. 

The  results  for  the  forward  velocity  sensor  failure  are  representative  only  for  the  sensors  that 
have  very  weak  cross-coupling  with  other  sensors:  specifically,  the  angle  of  attack  and  sideslip  angle 
sensors.  The  other  sensors  have  stronger  cross-coupling.  This  cross-coupling  is  not  through  the 
output  matrix,  but  through  the  state  transition  matrix  (i.e.,  the  system  dynamics  and  kinematics)  and 
Kalman  filter  gain  matrix.  If  the  Kalman  filter  assumes  that  the  pitch  rate  measurement  is  accurate 
(which  the  fully  functional  Kalman  filter  does),  the  state  estimates  are  updated  using  this 
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Table  6 


Kalman  Filter  Residual  Statistical  Comparison  for  the  Case  of  a  Pitch  Rate  Sensor  Failure. 


Residual 

Element 

Computed 

Average 

of^ 

Ensemble 

Temporal 

Average 

of? 

Temporal 

Ensemble 

Average 

of? 

Computed 
o  of  ? 

Ensemble 
Temporal 
o  of  ? 

Temporal 
Ensemble 
o  of  ? 

Forward  Velocity 

0 

7.25  E-2 

2.29  E-1 

3.003 

2.98 

2.97 

Angle  of  Attack 

0 

3.12  E-4 

-6.13  E-4 

2.96  E-2 

2.95  E-2 

2.92  E-2 

Pitch  Rate 

0 

1.56  E-2 

1.04  E-2 

3.46  E-2 

3.47  E-2 

3.39  E-2 

Pitch  Angle 

0 

-1.08  E-3 

-6.76  E-3 

2.92  E-3 

2.94  E-3 

2.94  E-3 

Sideslip  Angle 

0 

7.02  E-4 

-3.04  E-4 

1.44  E-1 

1.43  E-1 

1.43  E-1 

Roll  Rate 

0 

2.67  E-4 

3.57  E-4 

3.46  E-2 

3.43  E-2 

3.46  E-2 

Roll  Angle 

0 

-3.49  E-5 

-2.79  E-4 

3.94  E-3 

3.96  E-3 

3.87  E-3 

Yaw  Rate 

0 

1.37  E-4 

-4.86  E-3 

5.32  E-2 

5.38  E-2 

5.31  E-2 

Yaw  Angle 

0 

1.20  E-5 

-4.48  E-5 

3.22  E-3 

3.20  E-3 

3.16  E-3 

measurement,  even  if  they  are  wrong!  Thus,  the  other  states  are  a)rruptcd  by  the  incorrect 
measurement  of  this  single  element  of  the  measurement  vector.  The  residuals  not  only  reflect  the 
difference  between  the  actual  measurement  and  the  failed  measurement,  but  also  differences  will 
occur  between  the  coupled  measurements  and  the  filter  prediction  of  those  measurements.  This  is 
clearly  evident  in  Figure  18,  where  the  pitch  angle  element  of  the  residual  is  impacted  as  well  as  the 
pitch  rate  element.  The  average  and  standard  deviation  of  the  difference  between  the  computed 
residual  mean  and  the  actual  residual  are  tabulated  in  Table  6. 
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Figure  18.  Residual  and  Computed  Residual  Mean  and  Standard  Deviation  of  the 
Fully  Functional  Kalman  Filter  with  a  Mismodeled  Pitch  Rate  Sensor. 
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4.3  Single  Residual  Kalman  Filter  Bank 


In  Section  3.2.5  we  developed  the  framework  for  an  MMAE  structure  that  uses  the  residual 
vector  from  a  single  Kalman  filter  with  a  linear  transform  to  produce  the  equivalent  residual  from 
another  Kalman  filter.  In  this  section  we  present  a  comparison  of  the  residual  from  a  fully 
implemented  Kalman  filter  with  the  computed  equivalent  of  that  residual.  The  application  that  we 
studied  had  an  MMAE  with  16  Kalman  filters.  Extensive  simulation  was  conducted,  since  each  of 
the  16  filters  could  be  used  to  compute  the  residuals  from  any  of  the  other  16  filters  for  any  of  the  16 
possible  failure  scenarios.  This  results  in  4096  test  cases  (16^),  so  we  have  elected  to  present  a 
representative  sample  of  three  test  cases. 

4.3.1  No  Mismodeling  rNo  Failure).  First,  we  present  the  test  case  in  which  the 

true  system  has  no  flight  conu-ol  failure  and  we  compare  the  actual  residual  from  the  Kalman  filter 
that  assumed  a  fully  functional  aircraft  with  its  equivalent  residual  computed  using  the  residual  from 
the  Kalman  filter  that  assumed  a  left  elevator  failure.  This  is  a  somewhat  arbitrary  choice  to 
demonstrate  the  equivalency  between  the  usual  filter  bank  structure  of  several  parallel  Kalman  filters, 
and  a  filter  "bank"  structure  that  uses  a  single  Kalman  filter  with  several  linear  transforms  to  produce 
equivalent  residuals.  We  expect  that  the  usual  implementation  of  this  equivalent  structure  would 
consist  of  the  fully  functional  model  based  Kalman  filter  with  the  equivalent  residuals  computed  for 
the  Kalman  filters  that  assume  the  various  failures.  Figure  19  shows  what  appears  to  be  a  set  of 
single  plots  of  the  fully  functional  Kalman  filter  residual  elements,  but  these  plots  are  actually  the 
equivalent  residual  plotted  over  the  actual  residual.  It  appears  to  be  a  single  plot  because  the  two 
residuals  are  equivalent.  We  computed  the  difference  between  the  actual  residual  and  the  equivalent 
residual,  and  temporally  averaged  this  difference.  The  results  are  tabulated  in  Table  7.  Note  that  for 
some  of  the  residual  elements,  the  temporal  average  of  the  difference  between  the  actual  residual  and 
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Figure  19.  The  Fully  Functional  Kalman  Filter  Residual  and  the  Equivalent  Residual 
Computed  Using  the  Left  Elevator  Kalman  Filter  Residual  for  a  No-Failure  Test  Case. 
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Table  7 


Temporal  Average  of  the  Difference  Between  the  Fully  Functional  Kalman  Filter  Residual 
and  the  Equivalent  Residual  Computed  Using  the  Left  Elevator  Failure  Filter  Residual 

for  a  No-Failure  Test  Case. 


Residual 

Temporal  Average  of 

Element 

Difference 

Forward  Velocity 

-1.376  E-15 

Angle  of  Attack 

4.342  E-19 

Pitch  Rate 

3.492  E-17 

Pitch  Angle 

-5.926  E-18 

Sideslip  Angle 

<  10-^* 

Roll  Rate 

o 

V 

Roll  Angle 

<  10-^' 

Yaw  Rate 

<  io-2> 

Yaw  Angle 

<  10-2' 

equivalent  residual  is  less  than  the  precision  of  the  simulation  software  (10'^').  and  in  all  cases  the 
difference  is  many  orders  of  magnitude  less  than  the  two  individual  residuals.  Obviously,  the  actual 
residual  and  the  computed  equivalent  residual  are  equivalent. 

In  Section  3.2.5,  we  proposed  using  these  equivalent  residuals  to  construct  the  Standard 
Kalman  Filter  Bank,  where  we  are  using  the  fully  functional  model  based  Kalman  filter  as  the  source 
filter  for  this  implementation.  We  verified  this  construct  in  Figure  20  where  we  show  the  residual 
from  the  elevator  failure  Kalman  filter  along  with  its  equivalent  residual  that  was  computed  using  the 
fully  functional  Kalman  filter.  This  test  case  matches  the  structure  shown  in  Figure  6  of  Section 
3.2.5,  for  the  case  of  an  elevator  failure.  The  tabulated  differences  between  the  actual  residual  and 
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Figure  20.  The  Left  Elevator  Kalman  Filter  Residual  and  the  Equivalent  Residual 
Computed  Using  the  Fully  Functional  Kalman  Filter  Residual  for  a  Left  Elevator  Failure 

Test  Case. 
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Table  8 


Temporal  Average  of  the  Difference  Between  the  Left  Elevator  Failure  Kalman  Filter  Residual 
and  the  Equivalent  Residual  Computed  Using  the  Fully  Functional  Kalman  Filter  Residual 

for  a  Left  Elevator  Failure  Test  Case. 


Residual 

Temporal  Average  of 

Element 

Difference 

Forward  Velocity 

1.372  E-15 

Angle  of  Attack 

-4.516  E-19 

Pitch  Rate 

-3.503  E-17 

Pitch  Angle 

5.912E-18 

Sideslip  Angle 

<  10-^' 

Roll  Rate 

<  io-2‘ 

Roll  Angle 

<  10-^' 

Yaw  Rate 

<  io-^‘ 

Yaw  Angle 

<  10-^' 

the  equivalent  residual  are  in  Table  8;  again,  the  equivalence  of  the  two  forms  of  the  residual  is 
evident. 

4.3.2  Mismodeled  Control  Input  Matrix  (Actuator  Failure).  Next,  we 

present  the  test  case  in  which  the  true  system  has  an  actuator  failure,  specifically  the  left  elevator,  and 
we  compare  the  actual  residual  from  the  Kalman  filter  that  uses  a  fully  functional  aircraft  model,  with 
its  equivalent  residual  computed  using  the  residual  from  the  Kalman  filter  that  assumed  a  left 
elevator  failure.  Figure  21  again  shows  what  appears  to  be  a  set  of  single  plots  of  the  fully  functional 
Kalman  filter  residual,  but  these  plots  are  actually  the  equivalent  residual  plotted  over  the  actual 
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Computed  Using  the  Left  Elevator  Failure  Kalman  Filter  Residual 


for  a  Left  Elevator  Failure  Test  Case. 
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1. 


residual.  We  also  tabulated  the  difference  between  the  actual  residual  and  the  equivalent  residual  in 
Table  9.  Once  again,  these  residuals  are  equivalent. 

To  verify  the  structure  developed  in  Section  3.2.5,  we  show  the  actual  residual  from  the 
Kalman  filter  that  assumes  a  left  aileron  failure  when  the  actual  failure  is  a  left  elevator,  along  with 
its  equivalent  residual  computed  using  the  fully  functional  Kalman  filter.  This  comparison  is  made  in 
Figure  22,  and  the  differences  are  tabulated  in  Table  10.  Clearly,  the  equivalent  residual  and  the 
actual  residual  are  identical. 


Table  9 

Temporal  Average  of  the  Difference  Between  the  Fully  Functional  Kalman  Filter  Residual 
and  the  Equivalent  Residual  Computed  Using  the  Left  Elevator  Failure  Kalman  Filter  Residual 

for  a  Left  Elevator  Failure  Test  Case. 


Residual 

Temporal  Average  of 

Element 

Difference 

Forward  Velocity 

1.372  E-15 

Angle  of  Attack 

-4.5 16  E- 19 

Pitch  Rate 

-3.503  E-17 

Pitch  Angle 

5.912  E-1 8 

Sideslip  Angle 

<  10" 

Roll  Rate 

<  10" 

Roll  Angle 

<  10" 

Yaw  Rate 

<  10" 

Yaw  Angle 

<  10" 
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Computed  Using  the  Fully  Functional  Kalman  Filter  Residual 
for  a  Left  Elevator  Failure  Test  Case. 
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Table  10 


Temporal  Average  of  the  Difference  Between  the  Left  Aileron  Kalman  Filter  Residual 
and  the  Equivalent  Residual  Computed  Using  the  Fully  Functional  Kalman  Filter  Residual 

for  a  Left  Elevator  Failure  Test  Case. 


Residual 

Temporal  Average  of 

Element 

Difference 

Forward  Velocity 

<  io-^‘ 

Angle  of  Attack 

<  10-^' 

Pitch  Rate 

<  10-^' 

Pitch  Angle 

<  io-'‘ 

Sideslip  Angle 

-1.255  E-18 

Roll  Rate 

-6.940  E-18 

Roll  Angle 

-5.321  E-17 

Yaw  Rate 

-1.931  E-17 

Yaw  Angle 

-3.987  E-17 

4.3.3  Mismodeled  Output  Matrix  (Sensor  Failure).  Finally,  we  examine  the 

test  case  in  which  the  true  system  has  a  sensor  failure,  specifically  the  pitch  angle  sensor,  and 
compare  the  actual  residual  from  the  Kalman  filter  that  assumed  a  pitch  rate  sensor  failure  with  its 
equivalent  residual  computed  using  the  residual  from  the  Kalman  filter  that  assumed  a  fully 
functional  aircraft.  Note  that  the  actual  failure,  the  assumed  failure  for  the  output  filter  residual,  and 
the  chosen  source  Filter  are  nonmatching  (three-way  nonmatching)  to  help  demonstrate  the  level  of 
viability  of  this  approach.  We  have  chosen  the  source  filter  to  be  the  fully  functional  filter  to  verify 
the  structure  shown  in  Figure  6  of  Section  3.2.5.  As  in  the  previous  cases,  Figure  23  shows  what 
appears  to  be  a  set  of  single  plots  of  the  fully  functional  Kalman  filter  residual,  but  these  plots  are 
actually  the  equivalent  residual  plotted  over  the  actual  residual.  Once  again,  it  appears  to  be  a  single 
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Figure  23.  The  Pitch  Rate  Failure  Kalman  Filter  Residual  and  the  Equivalent  Residual 
Computed  Using  the  Fully  Functional  Kalman  Filter  Residual 
for  a  Pitch  Angle  Sensor  Failure  Test  Case. 


Table  11 

Temporal  Average  of  the  Difference  Between  the  Pitch  Rate  Failure  Kalman  Filter  Residual 
and  the  Equivalent  Residual  Computed  Using  the  Fully  Functional  Kalman  Filter  Residual 

for  a  Pitch  Angle  Failure  Test  Case. 


Residual 

Temporal  Average  of 

Element 

Difference 

Forward  Velocity 

-2.194  E-15 

Angle  of  Attack 

3.323  E-19 

Pitch  Rate 

-3.686  E-19 

Pitch  Angle 

4.966  E-18 

Sideslip  Angle 

<  10'' 

Roll  Rate 

<  10" 

Roll  Angle 

1.169  E-19 

Yaw  Rate 

-2.575  E-21 

Yaw  Angle 

-3.287  E-21 

plot  because  the  two  residual  are  equivalent.  The  differences  are  tabulated  in  Table  1 1 ,  and  again 
demonstrate  the  equivalency. 

4.3.4  Overall  Equivalency.  The  average  difference  for  all  4096  test  cases  was 

computed  and  tabulated  in  Table  12.  These  values  indicate  that  the  equivalent  residuals  are  identical 
(within  round-off  error)  to  the  actual  Kalman  filter  residuals. 

4.3.5  Generalized  Likelihood  Ratio  Class  Equivalence.  We  have  [39] 

noticed  the  striking  similarity  between  the  filter  "bank"  structure  of  Figure  6  and  Generalized 
Likelihood  Ratio  (GLR)  based  failure  detection  algorithms.  The  GLR  structure  [66]  uses  a  Kalman 
filter  residual  as  the  input  to  a  set  of  matched  filters,  and  then  the  output  of  these  matched  filters  are 
used  by  a  hypothesis  testing  algorithm  to  determine  the  failure  status.  The  matched  filters  are  based 
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Table  12 


Overall  Temporal  Average  of  the  Difference  Between  an  Actual  Residual 
and  the  Equivalent  Residual,  Averaged  over  all  4096  Test  Cases. 


Residual 

Temporal  Average  of 

Element 

Difference 

Forward  Velocity 

9.633  E-17 

Angle  of  Attack 

-1.378  E-18 

Pitch  Rate 

-1.518  E-17 

Pitch  Angle 

7.768  E-18 

Sideslip  Angle 

1.172  E-18 

Roll  Rate 

2.686  E-19 

Roll  Angle 

2.730  E-18 

Yaw  Rate 

-3.139  E-18 

Yaw  Angle 

1.576  E-18 

on  the  designer's  expectation  of  the  time  history  of  the  residual  for  certain  failures.  For  example,  the 
set  of  filters  would  model  a  no-failure  hypothesis,  a  jump  in  the  residual  mean  (representing  a  sudden 
additive  bias),  and  a  ramp  in  the  residual  mean  (representing  a  slow  changing  additive  bias). 

The  results  from  this  section  have  shown  that  the  linear  transforms  of  Figure  6  are 
equivalent  to  implementing  a  complete  Kalman  filter.  Another  interpretation  of  this  result  is  that  the 
linear  transforms  operating  on  a  single  Kalman  filer  output  can  be  viewed  as  the  best  possible 
"matched  filters"  for  the  GLR  failure  identification  method.  These  linear  transforms  yield  the  actual 
impact  of  the  hypothesized  failure  on  the  Kalman  filter  residual,  versus  the  designer's  ad  hoc 
representation  of  what  that  impact  might  be.  This  interpretation  establishes  an  equivalence  class 
relationship  between  MMAE  failure  detection  and  isolation  systems  and  the  best-possible  GLR 
failure  detection  and  isolation  systems.  This  equivalence  class  relationship  needs  to  be  established 
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using  a  more  structured  theoretical  approach.  This  appears  as  one  of  the  recommendations  for  future 
research,  in  Section  5.6. 

4.4  Failure  Detection  Performance 

We  present  the  failure  identification  performance  for  the  three  MMAE  structures  that  were 
developed  in  Chapter  3.  The  first  structure,  shown  in  Figure  24,  is  the  Standard  Kalman  Filter  Bank 
(SKFB)  with  a  Standard  Hypothesis  Testing  Algorithm  (SHTA).  The  Kalman  filter  equations  for  the 
SKFB  were  presented  in  Section  3.3.1,  and  the  SHTA  algorithm  was  presented  in  Section  3.3.1.  An 
alternative  structure,  shown  in  Figure  25,  was  developed  in  Section  3.2.5  and  was  shown  to  be 


equivalent  to  the  SKFB  in  Section  4.3.  Since  the  structures  are  equivalent,  the  failure  identification 


Figure  24.  Multiple  Model  Adaptive  Estimation  Algorithm 
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Figure  25.  Alternative  Multiple  Model  Adaptive  Estimation  Filter  Bank  Structure 

Using  Equivalent  Residuals 


performance  of  the  first  structure  will  be  identical  to  the  failure  identification  performance  of  the 
equivalent  sUaicture.  Therefore,  we  did  not  test  the  failure  identification  performance  of  this 
equivalent  structure.  The  second  structure,  shown  in  Figure  26,  is  the  Residual  Correlation  Kalman 
Filter  Bank  (RCKFB),  developed  in  Section  3.2.6,  with  a  SHTA.  The  last  structure  is  the  SKFB  with 
a  Neyman-Pearson  Hypothesis  Testing  Algorithm  (NPHTA)  that  was  developed  in  Section  3.3.2.  In 
Section  3. 3.2.1,  we  showed  that  performing  a  hypothesis  test,  using  a  likelihood  ratio,  on  multiple 
residuals  was  equivalent  to  performing  multiple  hypothesis  tests  on  a  single  residual.  Thus,  we  can 
replace  the  K  filters  in  the  Kalman  filter  bank,  with  a  single  filter  and  perform  a  Neyman-Pearson 
Hypothesis  Test  on  the  residual  from  this  single  filter.  If  the  MMAE  weighted  state  estimate  is 
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Figure  26.  Multiple  Model  Adaptive  Estimation  Algorithm  using  a 
Residual  Correlation  Kalman  Filter  Bank 


needed,  we  would  need  to  implement  either  the  full  SKFB,  or  the  equivalent  filter  bank  so  that  we 
have  all  of  the  state  estimates  from  the  various  Kalman  filters.  Thus,  the  SKFB-NPHTA  structure  is 
equivalent  to  the  structure  shown  in  Figure  27.  Since  we  do  not  need  the  state  estimates  to  test  the 
failure  identification  performance  of  this  structure,  we  did  not  implement  the  linear  transforms.  We 
chose  to  construct  the  SKFB-NPHTA  sumcture  using  the  fully  functional  aircraft  model  based 
Kalman  filter.  The  failure  identification  performance  of  the  SKFB-SHTA  is  presented  in  Section 
4.4.1,  the  RCKFB-SHTA  performance  in  Section  4.4.2,  and  the  performance  for  the  SKFB-NPHTA 
is  shown  in  Section  4.4.3.  Section  4.4.4  is  a  comparison  of  the  failure  identification  performances 
for  these  three  structures.  The  last  possible  combination  of  these  structures,  namely  the  RCKFB- 
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Figure  21.  Equivalent  Standard  Kalman  Filter  Bank  (SKFB)  with  a 
Neyman-Pearson  Hypothesis  Testing  Algorithm  Structure 


-NPHTA,  is  a  viable  (and  promising)  algorithm,  but  it  was  not  evaluated  because  of  lime  limitations. 

In  Section  3.2.4  we  found  that  the  mean  of  the  residual  will  have  elements  of  the  input  for 
the  case  of  an  actuator  failure,  and  for  a  sensor  failure  the  residual  will  have  elements  of  the  state 
estimates.  If  we  use  a  sinusoid  for  the  input,  then  the  sinusoid  input  elements  cause  the  residual  to  be 
a  sinusoid.  This  sinusoid  appears  as  a  spike  in  the  power  spectral  density  at  the  frequency  of  the 
input,  which  was  shown  in  Figure  9.  Therefore,  since  we  can  stipulate  the  frequency  of  the  system 
input  (since  we  can  generate  the  dither  input  ourselves),  we  can  use  the  spectral  content  at  that 
frequency  to  indicate  the  presence  of  an  actuator  failure,  which  was  the  basis  for  developing  the 
RCKFB.  However,  we  cannot  readily  stipulate  the  frequency  of  the  state  estimates,  unless  extensive 
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simulation  is  done  to  find  the  particular  input  that  causes  sinusoidal  state  estimates.  Since  the  state 
estimates  are  not  necessarily  sinusoidal,  they  may  not  be  well  correlated  and  we  might  need  to 
estimate  the  power  spectral  density  across  the  entire  spectrum  to  try  to  find  indications  of  the 
presence  of  a  sensor  failure.  Figure  1 8  shows  a  residual  from  a  filter  with  mismodeled  output  matrix. 
Clearly,  the  components  of  this  residual  are  not  sinusoidal.  We  believe  that  this  technique  would  not 
work  well  for  sensor  failures  unless  the  appropriate  input  is  used  that  causes  the  states  estimates  to 
become  strongly  correlated  in  a  readily  distinguishable  manner. 

For  this  reason,  we  implemented  the  RCKFB-SHTA  structure  using  only  actuator  failure 
models.  The  failure  identification  of  the  other  structures  was  initially  evaluated  for  all  the 
hypothesized  failures  (including  sensor  failures),  and  then  the  set  of  hypotheses  was  reduced  to  just 
the  actuator  failures  so  we  could  compare  the  failure  identification  performance  between  the  various 
structures.  This  reduced  set  of  hypotheses  significantly  reduced  the  computations  required  to  analyze 
the  failure  identification  performance,  so  we  were  able  to  compare  the  failure  identification 
performance  of  the  various  structures  at  several  different  system  input  levels.  We  also  restricted  the 
SKFB-SHTA  and  SKFB-NPHTA  (to  be  presented  in  Section  4.4.1)  structures  to  actuator  failure 
hypotheses  only,  to  make  a  fair  comparison  of  failure  identification  algorithms.  The  usual 
implementation  of  these  structures  includes  sensor  failure  hypotheses. 

Thus,  in  the  following  sections  you  will  first  find  the  failure  identification  performance  of 
each  of  the  structures,  and  then  a  comparative  evaluation  of  their  performances  for  the  different  types 
of  actuator  failures.  In  their  individual  failure  identification  performance  sections,  we  present  the 
overall  failure  identification  performance  of  the  SKFB-SHTA  and  the  SKFB-NPHTA,  including 
sensor  failures,  but  only  the  actuator  failure  identification  for  the  RCKFB-SHTA.  The  actuator 
failure  identification  performance  as  a  function  of  input  level  is  presented  for  all  three  structures,  and 
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then  the  presentation  of  the  actuator  performance  is  rearranged  to  facilitate  the  comparison  of  the 
performance  of  the  various  structures. 

4.4.1  SKFB-SHTA  Failure  Identification  Performance.  The  failure 

identification  performance  of  the  SKFB-SHTA  structure  for  this  particular  application  was 
previously  researched  [19, 41]  and  some  optimization  of  the  modifications  presented  in  Section  3.3.1 
was  accomplished.  We  present  a  summary  of  those  results  to  facilitate  performance  comparisons  to 
the  other  structures. 

The  failure  identification  performance  for  a  particular  failure  (a  left  elevator  failure)  is 
presented  in  Figure  28.  Note  that  the  conditional  probability  for  the  fully  functional  hypothesis 
essentially  mirrors  the  conditional  probability  for  the  actual  failure  (the  left  elevator  failure).  Similar 
results  were  observed  for  all  other  hypothesized  failures,  except  for  two.  The  failure  identification 
was  fairly  slow  for  rudder  failures  (the  left  one  is  shown  in  Figure  29)  and  the  forward  velocity 
sensor  failure  (shown  in  Figure  30).  We  present  these  failures  to  show  two  reasons  for  slow  failure 
identification  performance.  The  rudder  failure  identification  is  slow  because  several  other 
hypotheses  appear  to  be  equally  as  viable  as  the  rudder  failure  hypothesis.  Note  that  die  roll  angle 
sensor  failure  and  sideslip  angle  sensor  failures  appear  to  be  almost  as  probable  as  the  fully 
functional  hypothesis,  until  the  rudder  failure  hypothesis  is  identified.  Previous  research  [19]  has 
shown  that,  unless  the  filters  are  properly  tuned  as  we  described  in  Section  3. 3. 1.4,  this  structure  has 
a  bias  toward  occasionally  generating  a  false  alarm  for  sensor  failures,  which  is  precisely  what  is 
observed  in  this  case.  The  forward  velocity  sensor  failure  performance,  Figure  30,  shows  a  different 
situation.  This  sensor  has  much  stronger  measurement  noise  than  all  the  other  sensors,  so  the  effects 
of  this  particular  sensor  failure  have  to  be  quite  large  for  the  failure  to  be  identified.  Figure  17  shows 


139 


Standard  Hypothesis  Testing  Algorithm  (SHTA)  for  a  Left  Elevator  Failure. 
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Figure  29.  Failure  Identification  of  the  Standard  Kalman  Filter  Bank  (SKFB)  with  a 
Standard  Hypothesis  Testing  Algorithm  (SHTA)  for  a  Left  Rudder  Failure. 
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Figure  30.  Failure  Identification  of  the  Standard  Kalman  Filter  Bank  (SKFB)  with  a 
Standard  Hypothesis  Testing  Algorithm  (SHTA)  for  a  Forward  Velocity  Sensor  Failure. 
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that  it  takes  a  few  seconds  for  the  effects  of  the  sensor  failure  to  cause  the  residual  to  grow.  In  this 
case  the  noise  of  the  sensor  is  so  great  that  it  masks  the  residual  growth  until  it  is  quite  large. 

The  results  from  these  three  and  all  other  failure  cases  can  be  summarized  by  observing  only 
the  conditional  probability  of  the  particular  failure  that  is  actually  occurring,  since  the  fully  functional 
hypothesis  mirrors  this  probability,  and  all  other  hypothesized  failures  are  near  zero  throughout  the 
simulation.  The  overall  performance  of  the  SKFB-SHTA  failure  identification  performance  can  be 
summarized  by  showing  the  conditional  probabilities  for  each  of  the  failure  hypotheses  on  a  series  of 
plots,  shown  in  Figure  31,  where  each  plot  shows  the  SKFB-SHTA  failure  identification  performance 
for  a  particular  failure.  Note  that  the  left  elevator  failure  identification  performance  (Figure  28)  is 
summarized  in  the  left  elevator  plot  of  Figure  31,  the  left  rudder  failure  identification  performance 
(Figure  29)  is  shown  in  the  left  rudder  plot  of  Figure  3 1 ,  and  the  forward  velocity  sensor  failure 
identification  performance  (Figure  30)  is  shown  in  the  forward  velocity  plot  of  Figure  3 1 . 

One  simple  method  of  declaring  a  failure  status  is  to  choose  a  probability  threshold  that  is 
used  as  a  declaration  trigger.  The  conditional  probabilities  are  compared  to  this  threshold,  and  when 
the  threshold  is  exceeded,  the  hypothesis  corresponding  to  that  conditional  probability  is  declared  as 
the  failure  status.  For  example,  if  the  threshold  was  set  to  a  probability  of  0.5  and  a  left  elevator 
failure  occurred,  when  the  left  elevator  conditional  probability  exceeds  0.5,  then  a  failed  left  elevator 
would  be  declared.  Figure  3 1  shows  that  this  would  occur  at  about  1 .5  seconds  into  the  simulation, 
or  about  0.5  seconds  after  the  failure  actually  occurred.  Using  this  threshold,  we  averaged  the 
actuator  failure  performance  of  this  MMAE  structure  over  5  Monte  Carlo  runs,  and  then  averaged 
the  left  and  right  actuator  failure  performance  results  to  get  an  average  over  10  test  samples.  We 
averaged  the  left  and  right  actuator  failure  performances  because  we  assume  that  they  are 
independent.  This  assumption  is  based  on  the  orthogonal  system  input  (Figure  13),  where  the  left 
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Figure  31.  Overall  Failure  Identification  Performance  of  a  Standard  Kalman  Filter  Bank 
(SKFB)  with  a  Standard  Hypothesis  Testing  Algorithm  (SHTA). 
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actuator  input  is  almost  exactly  out-of-phase  with  the  right  input  (a  slight  offset  from  exactly  out-of- 
phase  was  needed  to  allow  some  excitation  of  the  sensors),  thus  the  left  and  right  inputs  are 
orthogonal.  The  averaging  was  done  for  a  various  input  strengths,  starting  with  the  one  shown  in 
Figure  13,  and  reducing  it  with  an  input  multiplier,  down  to  1/100  of  the  original  control  input.  The 
results  are  tabulated  in  Table  13,  where  the  decision  time  is  the  amount  of  time  from  the  occurrence 
of  the  failure  to  when  the  conditional  probability  of  the  hypothesis  crosses  the  probability  threshold 
of  0.5.  At  very  small  inputs,  the  conditional  probability  grew  so  slowly,  that  when  they  approached 
the  threshold,  the  randomness  of  the  probabilities  would  cause  it  to  cross  the  threshold  a  few  times. 
For  these  cases,  the  decision  time  is  the  amount  of  time  from  the  occurrence  of  the  failure  to  when  the 
conditional  probability  crosses  the  probability  threshold  for  the  last  time  and  stays  there.  These 
results  are  graphically  shown  in  Figure  32,  where  the  input  multiplier  is  plotted  on  a  logarithmic 
(base  10)  scale. 

Note  the  "  entries  in  Table  13  and  the  dashed  line  portion  of  the  plots  in  Figure  32.  These 
indicate  where  false  alarms  were  observed,  with  an  increase  in  the  number  of  false  alarms  as  the 
input  was  decreased.  The  dashes  are  to  show  where  the  performance  may  be  deemed  unacceptable 
due  to  the  increase  in  the  false  alarm  rate. 

Figure  32  shows  that  failure  performance  drops  off  dramatically  as  the  input  decreases.  This 
is  particularly  evident  for  the  elevator  and  aileron  failures,  where  identification  times  are  beyond  the 
termination  of  the  simulation  (7  seconds  after  the  failure)  for  input  multipliers  of  less  than  0.025. 
These  results  make  good  sense,  because  stronger  inputs  would  cause  larger  residuals  from  the 
incorrect  Kalman  filters  when  the  failure  occurs.  Likewise,  smaller  inputs  will  produce  relatively 
smaller  residuals  when  the  failure  occurs,  so  the  conditional  probability  of  the  failure  hypothesis 
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Table  13 


SKFB-SHTA  Actuator  Failure  Identification  Performance,  10-Test  Sample  Average. 


Input 

Elevator  Failures 

Aileron  Failures 

Rudder  Failures 

Multiplier 

(sec) 

(sec) 

(sec) 

1.00 

0.42 

0.31 

2.29 

0.25 

0.51 

0.41 

1.04 

0.1 

0.96 

0.98 

1.11 

0.075 

1.38 

2.07 

1.29 

0.0625 

2.07  - 

3.34  - 

1.23  - 

0.05 

3.28  - 

>7- 

2.495  - 

0.025 

>7- 

>7- 

4.98  - 

Figure  32.  Actuator  Failure  Identification  Performance  for  the  Standard 
Kalman  Filter  Bank  (SKFB)  with  a  Standard  Hypothesis  Testing  Algorithm 

(SHTA). 
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requires  more  time  to  grow.  In  particular,  note  that  the  rudder  failure  identification  is  still  viable  at 
an  input  multiplier  level  of  0.025,  but  drops  off  quickly  after  that. 

Note  that,  as  the  rudder  input  is  increased,  the  failure  identification  time  actually  increases, 
which  does  not  correspond  to  the  anticipated  results  of  better  identification  performance  for  stronger 
inputs.  We  showed  earlier  that  this  structure  had  trouble  identifying  rudder  failures,  possibly 
because  of  cross-coupling  between  axes  causing  misidentification  of  a  sensor  failure.  This  is  not  the 
case  here,  since,  at  this  point,  we  have  limited  the  SKFB  to  actuator  failure  hypotheses  to  make  a  fair 
comparison  with  the  RCBCFB  structure.  There  are  two  reasons  that  this  occurs,  one  is  the  stonger 
rudder  input  (Figure  13),  relative  to  the  elevator  and  aileron  inputs,  and  the  other  is  the  timing  of  the 
failure. 

To  illustrate  these  reasons,  we  have  plotted  the  likelihood  quotient,  from  Eq  (99),  from 
each  of  the  Kalman  filters  for  the  no-failure  case.  These  quotients  are  used  as  negative  exponents  in 
Eq  (103),  if  p  stripping  described  in  Section  3.3. 1.2  is  not  being  used,  and  Eq  (105),  if  p  stripping  is 
being  used  (we  implemented  the  SHTA  using  p  stripping  for  this  research),  and  are  key  factors  in 
computing  the  conditional  probabilities.  Large  quotients  produce  extremely  small  exponential 
factors  in  Eq  (103)  and  Eq  (105),  which,  when  multiplied  by  the  prior  conditional  probability,  tends 
to  decrease  the  elemental  conditional  probability  significantly.  These  elemental  conditional 
probabilities  are  normalized  by  the  sum  of  all  the  conditional  probabilities.  Good  failure 
identification  requires  that  the  conditional  probabilities  change  quickly  when  a  failure  occurs,  which 
in  turn  requires  that  the  quotient  for  the  correct  hypothesis  be  small  compared  to  the  other  quotients, 
so  that  its  conditional  probability  will  not  decrease  as  much  as  the  others.  This  is  clearly  indicated  in 
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Figure  33,  where  the  quotient  for  the  no-failure  Kalman  filter  is  much  smaller  than  the  other 
quotients. 

Using  the  plots  on  Figure  33,  we  see  that,  when  the  simulation  is  started,  the  strong  rudder 
input  causes  the  quotient  in  the  rudder  failure  Kalman  filter  to  grow  quite  large,  much  larger  than  the 
other  quotients.  When  the  right  rudder  failure  actually  occurs  at  one  second  into  the  simulation  in 
Figure  34,  the  rudder  quotient  is  nearly  at  its  maximum,  and  almost  five  times  larger  than  any  other 
quotient.  During  the  first  second  (i.e.,  with  no  failure),  the  large  rudder  input  causes  the  state 
estimates  to  diverge  rapidly  from  the  true  states,  thus  causing  a  large  residual,  which 
in  turn  causes  a  large  rudder  likelihood  quotient.  This  divergence  occurs  in  the  other  incorrect  filters 
as  well,  but  to  a  lesser  degree  due  to  the  smaller  inputs,  which  is  why  their  likelihood  quotients  are 
smaller.  When  the  failure  occurs,  the  state  estimates  in  the  rudder  failure  Kalman  filter  (now  the 
correct  filter)  take  longer  to  converge  to  the  true  states  than  a  filter  whose  state  estimates  were  closer 
to  the  true  states  when  it  became  the  correct  filter.  Thus  the  residual  will  take  longer  to  decrease,  for 
the  rudder  failure  case,  which  causes  the  likelihood  quotient  to  decrease  slowly,  as  shown  in  Figure 
34.  This  causes  the  conditional  probability  to  increase  slowly,  which  delays  the  failure  identification. 

The  timing  of  the  failure  is  also  important  because  of  the  interplay  between  the  various 
likelihood  quotients.  Figure  33  shows  that  at  one  second  into  the  simulation  the  rudder  likelihood 
quotient  is  at  its  maximum  and  the  other  likelihood  quotients  are  at  a  local  minimum.  Figure  34 
shows  that  this  causes  a  delay  in  the  deaease  of  the  rudder  likelihood  quotient,  and  the  increase  of 
the  other  likelihood  quotients,  which  delays  the  change  in  the  conditional  probabilities.  If  the  failure 
is  injected  at  4.2  seconds  into  the  simulation  instead,  the  rudder  likelihood  quotients  (Figure  35)  are 
approaching  their  local  minimum,  the  aileron  quotients  are  steady,  and  the  elevator  quotients  are 
starting  to  rise  to  a  local  maximum.  These  conditions  favor  a  faster  rudder  failure 
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Figure  33.  Kalman  Filter  Likelihood  Quotients  for  the  No-Failure  Case. 
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Figure  34.  Kalman  Filter  Likelihood  Quotients  for  the  Right  Rudder  Failure  Case. 
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Figure  35.  Kalman  Filter  Likelihood  Quotients  for  a  Right  Rudder  Failure 

Injected  at  4.2  Seconds. 
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identification  time  because  the  mdder  quotient  is  already  diminishing,  while  the  other  quotients  are 
increasing.  On  the  other  hand,  these  conditions  appear  to  be  unfavorable  if  an  aileron  failure  were  to 
occur  at  4.2  seconds,  since  we  found  that  the  aileron  failure  detection  time  was  much  greater,  under 
these  conditions,  than  we  found  for  the  failure  injection  times  of  1  and  2  seconds.  We  found  that  it 
took  about  2.4  seconds  to  identify  the  failure  clearly;  during  that  time  the  aileron  failure  conditional 
probability  and  the  fully  functional  conditional  probability  oscillated  between  each  other,  finally 
settling  onto  the  correct  hypothesis  at  about  2.4  seconds.  Clearly,  failure  timing  in  relation  to  the 
likelihood  quotients  can  strongly  influence  the  failure  identification  performance. 

We  conducted  a  the  same  type  of  analysis  that  was  used  to  generate  Figure  32,  but  we  did 
this  for  different  failure  injection  times.  The  failures  were  injected  at  1,  2,  and  4.2  seconds  into  the 
simulation.  The  input  multiplier  is  plotted  as  a  function  of  failure  identification  time  for  the  elevator 
(Figure  36),  aileron  (Figure  37),  and  rudder  (Figure  38).  The  shaded  area  represents  the  range  of 
identification  times  that  we  believe  could  occur  for  the  given  input  (shown  in  Figure  13),  but  with  the 
different  multipliers.  Figure  36  indicates  that  there  is  an  optimal  input  strength  for  which  the 
difference  between  the  maximum  (worst)  failure  identification  time  and  the  minimum  (best)  failure 
identification  time  is  minimized.  This  optimal  input  strength  seems  to  be  just  prior  to  the  dramatic 
dropoff  in  failure  identification  performance.  This  makes  good  sense  because,  as  the  input  decreases, 
the  relative  strengths  between  the  various  inputs  also  decrease,  so  the  relative  strengths  of  the 
likelihood  quotients  are  also  decreasing.  For  example,  let  the  inputs  be  20,  10,  and  15  degrees  for 
the  elevator,  aileron,  and  rudder,  respectively.  If  the  multiplier  is  .25,  then  the  inputs  become  5, 2.5, 
and  3.75  degrees.  If  the  multiplier  is  0.1,  then  the  inputs  are  2,  1,  and  1.5  degrees.  Obviously,  the 
range  between  the  inputs  is  decreasing,  which  would  cause  a  decrease  in  the  range  between  the 
likelihood  quotients. 
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Figure  36.  Elevator  Failure  Identification  Performance 
for  Various  Failure  Injection  Times. 


Figure  37.  Aileron  Failure  Indentification  Performance 
for  Various  Failure  Injection  Times. 
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Figure  38.  Rudder  Failure  Identification  Performance 
for  Various  Failure  Injection  Times. 


Figure  36  shows  excellent  elevator  failure  identification  performance  with  very  little 
difference  between  the  maximum  and  minimum  identification  times.  We  believe  this  is  due  not  only 
to  the  strength  of  the  elevator  likelihood  quotients,  but  also  the  oscillatory  behavior  that  is  clearly 
seen  in  Figure  33.  This  oscillatory  behavior  would  help  cause  the  likelihood  quotient  to  build  if  the 
hypothesis  is  incorrect,  and  diminish  quickly  if  the  hypothesis  is  correct.  These  observations  would 
argue  for  finding  a  "balanced"  input,  one  that  causes  the  oscillatory  behavior  seen  in  the  elevator 
quotients  in  Figure  33,  and  cause  equal  magnitude  likelihood  quotients  for  all  of  the  incorrect 
hypotheses,  in  order  to  get  good  detection  performance  for  ah  types  of  actuator  failures. 

Time  constraints  prevented  further  analysis  of  these  observations.  In  particular,  we  were 
not  able  to  compare  the  performance  of  other  structures  for  the  different  failure  injection  times,  but 
we  were  able  to  compare  the  other  structures'  performance  for  the  given  input,  shown  in  Figure  13. 
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The  comparison  of  the  failure  identification  performance  for  the  various  structures  will  be  presented 
in  Section  4,4.4. 

4.4.2  RCKFB-SHTA  Failure  Identification  Performance.  We  implemented 

the  Residual  Correlation  Kalman  Filter  Bank  (RCKFB)  with  the  Standard  Hypothesis  Testing 
Algorithm  (SHTA)  shown  in  Figure  26.  We  used  the  form  of  the  RCKFB,  shown  in  Eq  (72),  that 
exploits  the  fast  Fourier  transform,  but  ignores  the  cross-coupling  terms  of  the  residual.  These  cross¬ 
coupling  terms  are  not  small  in  many  cases.  For  example,  the  residual  for  a  left  elevator  failure 
(Figure  14)  shows  that  both  the  pitch  rate  and  pitch  angle  strongly  reflect  the  elevator  input.  Note 
that  the  pitch  rate  is  clearly  a  sinusoid,  and  would  be  strongly  correlated  with  itself.  The  pitch  angle 
is  also  a  sinusoid  at  tire  same  frequency,  so  the  cross-correlation  between  these  two  sinusoids  would 
be  quite  strong.  By  ignoring  these  cases  of  strong  cross-correlations,  we  are  implementing  an 
algorithm  that  does  not  exploit  any  of  the  information  available  from  these  cross  terms.  We 
considered  the  computational  cost  of  exploiting  these  cross  terms  to  be  unacceptable,  so  we 
implemented  the  algorithm  using  only  the  autocorrelation  of  the  individual  elements  of  the  residual. 

A  summary  of  the  failure  detection  performance,  following  the  format  described  in  the 
previous  section,  is  presented  in  Figure  39.  Time  limitations  prevented  investigating  the  sensor 
failure  detection  performance  of  this  structure. 

In  Section  3.2.6  we  noted  that  Kay  [24]  shows  that  this  algorithm  essentially  filters  the 
residual  with  a  bandpass  filter,  centered  at  the  frequency  of  interest  and  with  a  bandwidth  of  UN.  If 
the  bandwidth  of  this  filter  is  small,  the  other  input  frequencies  will  not  have  spectral  content  in  the 
bandwidth  of  this  bandpass  filter.  Thus,  we  need  to  collect  enough  data  samples  so  that  the 
bandwidth  of  the  filter  will  not  have  any  spectral  content  from  the  other  inputs.  The  more  data 
samples  that  we  need  to  collect,  the  longer  the  delay  will  be  from  the  time  when  the  sinusoid  is 
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Figure  39.  Actuator  Failure  Identification  Performance  of  a  Residual  Correlation  Kalman 
Filter  Bank  (RDKFB)  with  a  Standard  Hypothesis  Testing  Algorithm  (SHTA). 


present  (when  the  failure  occurs)  to  when  we  can  compute  the  spectral  estimate  of  this  filter  and 
detect  the  sinusoid.  Thus,  there  is  a  tradeoff  between  the  number  of  data  samples  needed  to 
distinguish  between  the  input  sinusoids  and  the  failure  detection  time.  This  would  dictate  a  wide 
spacing  of  the  input  sinusoids,  to  allow  the  bandpass  filters  to  have  larger  bandwidths,  so  fewer  data 
samples  are  needed.  Unfortunately,  most  mechanical  system  (like  aircraft)  have  very  narrow  low 
pass  system  response  bandwidths  (the  bandwidth  was  about  2.2  Hz  for  this  particular  aircraft).  The 
system  response  to  the  input  must  be  seen  in  the  residual  for  any  failure  identification  scheme  to 
work,  thus  we  need  to  specify  input  frequencies  that  are  within  the  system  response  bandwidth.  This 
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requires  that  the  inputs  be  located  close  together  in  the  frequency  spectrum,  which  requires  narrow 
bandpass  filters  for  this  structure,  which  in  turn  requires  more  data  samples  to  be  collected.  Our 
inputs  were  at  0.5, 1  and  2  Hz,  which,  when  normalized  by  dividing  by  the  sampling  frequency, 
become  0.01, 0.02,  and  0.04,  respectively.  The  two  closest  frequencies  were  0.5  Hz  away  from  each 
other,  or  0.01  in  normalized  frequency.  The  bandpass  filter  would  need  to  have  a  bandwidth  of  less 
than  1  Hz  (0.5  Hz  on  each  side  of  the  center  frequency),  which  is  0.02  in  normalized  frequency,  to  be 
able  to  distinguish  these  input  frequencies  from  each  other.  This  requires  100  data  samples  (2/iV  = 
0.02  =*  N=  100)  for  the  required  bandpass  filter  bandwidth,  which  is  the  value  of  N  chosen  for  this 
research.  This  value  caused  a  delay  in  the  detection  of  the  failure  because  at  least  half  of  the  100 
data  samples  (which  translates  to  1  second)  had  to  show  strong  correlation  at  the  input  sinusoid 
frequency  for  the  spectral  estimator  even  to  start  to  identify  the  presence  of  the  sinusoid.  Figure  39, 
when  compared  to  Figure  31,  shows  that  this  delay  caused  this  algorithm  to  take  about  one  second 
longer  to  identify  the  actuator  failures. 

We  noted  earlier  that  the  airaaft  model  assumed  complete  axis  decoupling,  which  can  help 
decrease  the  number  of  required  data  samples.  The  elevator  input  was  at  a  frequency  of  1  Hz,  while 
the  rudder  input  was  at  0.5  Hz  and  the  aileron  input  was  at  2  Hz.  Since  the  pitch  axis  would  not 
affect  the  roll  and  yaw  axes,  the  elevator  input  would  not  appear  in  the  spectra  for  the  longitudinal 
elements  of  the  residual  (roll  and  yaw).  Thus  the  elevator  input  frequency  did  not  need  to  be 
considered  when  we  were  designing  the  bandwidth  of  the  bandpass  filter.  Thus,  the  two  closest 
frequencies  were  actually  1.5  Hz  apart,  which  means  that  we  could  have  used  only  34  data  samples  to 
distinguish  between  the  rudder  and  aileron  input  frequencies.  Time  constraints  prevented  testing  the 
failure  identification  performance  using  other  values  of  N . 
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We  tested  the  failure  performance  of  this  structure  at  the  various  input  strengths  that  were 
used  for  the  SKFB-SHTA  testing.  These  results  are  tabulated  in  Table  14  and  graphically  presented 
in  Figure  40.  Note  that  the  results  are  an  average  of  only  two  Monte  Carlo  runs  with  both  left  and 
right  actuators,  for  a  total  of  four  test  samples.  We  observed  a  large  drop  in  failure  identification 
time  for  the  strong  rudder  inputs.  The  reason  for  this  phenomenon  is  illustrated  in  Figure  20.  Note, 
when  the  failure  occurs  at  one  second,  the  pitch  rate  and  pitch  angle  elements  of  the  residual  are  at 
their  peak  and  take  about  0.5  seconds  to  die  out.  During  this  decay,  these  waveforms  appear  to  be 
part  of  the  sinusoid  that  was  clearly  present  before  the  failure  occurred.  The  sinusoid  in  the  fully 
functional  Kalman  filter  (Figure  21)  is  also  delayed  by  about  0.5  second.  Thus,  the  left  elevator 
failure  Kalman  filter  residual  will  show  a  slowly  decaying  spectral  content  at  the  elevator  input 
frequency,  and  the  fully  functional  Kalman  filter  residual  will  show  a  slowly  building  spectral  content 
at  that  frequency.  Thus,  the  SHTA  will  not  identify  the  failure  until  the  decay  is  nearly  complete. 

The  important  feature  to  notice  in  Figure  40,  is  that  failure  identification  is  still 
accomplished  at  a  much  lower  input  level  than  accomplished  with  the  SKFB  algorithm,  as  seen  in 
Figure  32.  This  feature  will  be  highlighted  when  this  structure  is  compared  to  the  others  in  Section 
4.4.4. 

The  plots  of  the  likelihood  quotient.  Figure  33,  suggest  an  alternative  structure  to  the 
RCKFB.  We  note  that  the  likelihood  quotient  also  reflects  the  sinusoidal  nature  of  the  residual, 
especially  in  the  elevator  likelihood  quotients.  A  spectral  estimator  could  be  used  to  estimate  the 
power  spectral  density  of  these  likelihood  quotients,  rather  than  of  the  residuals.  Since  this  quotient 
is  a  scalar  that  is  formed  by  computing  the  scaled  squared  value  of  the  residual,  the  cross-coupling 
terms  would  be  included  in  this  structure.  Also,  the  spectral  estimator  would  be  fairly 
straightforward  to  implement,  due  to  the  decrease  in  the  dimension  of  the  input  to  the  estimator. 
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Table  14 


RCKFB-SHTA  Actuator  Failure  Identification  Performance,  4-Test  Sample  Average. 


Input 

Multiplier 

Elevator  Failures 
(sec) 

Aileron  Failures 
(sec) 

Rudder  Failures 
(sec) 

1.00 

2.55 

2.42 

6.35 

0.25 

3.10 

2.12 

4.52 

0.1 

1.85 

2.73 

3.15 

0.075 

1.72 

2.43 

2.26 

0.0625 

1.75 

2.26 

2.04 

0.05 

1.99 

2.29 

1.54 

0.025 

1.95 

2.08 

2.03 

0.016 

3.28  - 

2.81  - 

3.79  - 

0.01 

3.43  - 

4.43  - 

5.59  - 

Figure  40.  Actuator  Failure  Identification  Performance  for  the  Residual 
Correlation  Kalman  Filter  Bank  (RCKFB)  with  a  Standard  Hypothesis 
Testing  Algorithm  (SHTA). 
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Note,  that  by  using  this  likelihood  quotient  (the  scaled  squared  norm  of  the  residual)  we  are 
essentially  squaring  the  sinusoidal  components  of  the  residual,  which  doubles  the  frequency,  so  the 
likelihood  quotient  oscillates  at  twice  the  frequency  as  the  residual.  This  would  allow  us  to  use  fewer 
data  samples  since  the  bandwidth  of  the  bandpass  filters  could  be  twice  as  large  as  they  were  for  the 
RCKFB.  This  approach  shows  great  promise  for  improving  the  performance  of  residual  correlation 
based  Kalman  filter  banks. 

4.4.3  SKFB-NPHTA  Failure  Identification  Performance.  Using  the 

development  in  Section  3.3.2,  we  implemented  a  Neyman-Pearson  Hypothesis  Testing  Algorithm 
(NPHTA),  similar  to  the  structure  shown  in  Figure  27.  We  were  only  testing  the  failure  identification 
performance  of  the  structure,  so  we  did  not  implement  the  linear  transforms  that  are  needed  to 
produce  the  blended  state  estimate.  Thus  our  structure  simply  consisted  of  a  single  Kalman  filter  (the 
ftilly  functional  Kalman  filter)  with  a  NPHTA.  The  output  of  this  algorithm  is  a  declaration  of  the 
failure  status,  not  a  set  of  conditional  probabilities.  We  chose  to  compare  the  declared  failure  status 
to  the  true  failure  status  to  evaluate  the  performance  of  this  algorithm.  In  Figure  41  we  present  the 
failure  identification  performance  of  this  structure,  by  plotting  the  agreement  (denoted  by  1)  and 
disagreement  (denoted  by  0)  of  the  declared  failure  status  with  the  true  failure  status. 

The  design  parameters  for  the  NPHTA  are  the  probability  of  detection  and  probability  of 
false  alarm,  which  were  0.999  and  0.01  respectively,  for  this  research.  Performance  sensitivity  to  the 
choice  of  these  two  parameters  could  be  investigated  for  a  given  application,  but  these  are  chosen  as 
reasonable,  representative  values.  First  note  the  one  misidentification  in  the  right  aileron  plot  of 
Figure  4 1 .  The  fully  functional  hypothesis  was  momentarily  (for  about  0.02  sec)  declared  as  the 
failure  status,  when  in  fact  the  right  aileron  failure  was  the  uue  hypothesis.  We  attribute  this  to  the 
probability  of  false  alarm  setting.  Note  that  this  is  a  false  alarm  (using  the  definitions  of  Section  3.3) 
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Figure  41.  Overall  Failure  Identification  Performance  of  a  Standard  Kalman  Filter  Bank 
(SKFB)  with  a  Neyman-Pearson  Hypothesis  Testing  Algorithm  (NPHTA). 
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rather  than  a  missed  alarm,  since  the  primary  hypothesis  before  the  test  was  the  failed  state  but  the 
identification  algorithm  chose  the  unfailed  state.  From  the  time  of  failure  (1  sec)  to  the  end  of  the 
simulation  there  are  48  tests  of  the  right  aileron  failure  hypothesis  vs.  the  fully  functional  hypothesis. 
The  probability  of  false  alarm  setting  that  was  used  would  roughly  give  us  one  false  alarm  for  every 
100  tests,  so  to  get  one  false  alarm  for  a  set  of  48  tests  does  seem  reasonable. 

Sensor  failures  showed  mixed  results.  The  failure  of  sensors  that  have  very  little 
measurement  noise  was  identified  quite  readily,  but  the  two  sensors  that  have  the  most  measurement 
noise  show  several  false  alarms.  This  algorithm  was  developed  under  the  assumption  that  the 
covariance  of  the  test  statistic  was  the  same  for  all  of  the  hypotheses.  This  is  true  for  the  actuator 
failures,  but  not  necessarily  true  for  the  sensor  failures.  If  the  sensors  have  very  little  measurement 
noise,  then  the  covariance  of  the  test  statistic,  for  that  particular  hypothesis,  is  essentially  equal  to 
the  covariance  for  an  actuator  failure  hypothesis.  If  the  sensor  has  a  large  amount  of  measurement 
noise,  then  the  covariance  matrix  will  differ  greatly  from  the  covariance  matrix  for  an  actuator  failure 
hypothesis.  We  believe  that  the  different  covariances  are  the  reason  for  the  large  number  of  false 
alarms  for  these  sensor  failure  hypotheses. 

The  failure  identification  performance  for  this  structure  was  tested  using  the  range  of  inputs 
from  the  previous  sections,  and  the  results  are  tabulated  in  the  Table  15.  The  plot  of  those  results 
(Figure  42)  shows  that  this  structure  performs  quite  well,  especially  when  compared  with  the  SKFB- 
SHTA,  Figure  32,  and  RCKFB-SHTA,  Figure  40  (except  for  small  input  strengths).  Note  that  this 
algorithm  identifies  the  failure  much  more  quickly  than  the  other  structures,  until  the  input  reaches 
0.05  times  the  original  input,  and  then  the  performance  drops  off  quite  rapidly.  We  found  an 
increase  in  false  alarms  for  input  levels  below  the  0.05  multiplier.  It  appears  that  the  0.05  multiplier 
level  is  the  minimum  input  strength  needed  to  provide  good  failure  identification  performance.  Note 


162 


that  this  structure's  performance  remains  relatively  constant  with  input  strength,  until  it  reaches  this 
minimum  input. 

In  particular,  note  that  the  rudder  failure  identification  performance  does  mil  reflect  the  drop 
in  performance  as  the  input  strengthens,  tike  we  saw  for  the  other  two  structures.  The  basis  of  the 
SHTA  is  that  the  correct  residual  is  small  compared  to  the  residual  from  the  incorrect  filters.  We've 
seen  that,  if  the  failure  occurs  when  the  appropriate  residual  is  near  a  peak,  the  failure  identification 
is  delayed.  This  algorithm  computes  a  time  history  of  the  residual,  based  on  a  certain  hypothesis,  and 
then  compares  the  measurement  history  of  the  residual  to  the  computed  time  history  (if  the 
distinguishability  is  great  enough)  to  determine  which  hypothesis  is  correct.  Thus,  the  algorithm 
does  not  have  to  wait  until  the  residual  decays  before  it  can  identify  the  failure.  It  is  possible  that  this 
algorithm  could  suffer  from  the  failure  injection  timing  problems  that  were  described  for  the  SHTA. 
If  the  failure  occurs  Just  after  the  discrimination  measure  was  reset,  when  the  discrimination  measure 
grows  large  enough  to  trigger  a  test,  the  test  statistic  will  not  have  grown  as  large  as  it  would  have  if 
the  failure  had  actually  occurred  right  when  the  discrimination  measure  was  reset.  Thus,  the  NPHTA 
could  miss  the  change  in  hypothesis  (missed  detection)  when  the  test  is  conducted.  The 
discrimination  measure  and  test  statistic  would  be  reset,  and  they  would  build  appropriately  to 
identify  the  failure  on  the  next  test.  We  observed  that  actuator  failure  hypothesis  tests  were  triggered 
about  every  5-10  time  samples  and  sensor  failure  hypothesis  tests  were  occurring  about  every  1-5 
data  samples.  Thus,  it  is  possible  that  the  failure  identification  could  be  delayed  about  0.2  seconds 
beyond  the  best  detection  time.  This  compares  quite  favorably  with  the  range  in  detection  times  that 
were  found  using  the  SKFB-SHTA  structure,  shown  in  Figures  36,  37,  and  38.  A  similar  analysis  of 
the  timing  differences  needs  to  be  accomplished  for  this  structure,  in  the  future. 
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Table  15 


SKFB-NPHTA  Actuator  Failure  Identification  Performance,  10-Test  Sample  Average. 


Elevator  Failures 

Aileron  Failures 

Rudder  Failures 

(sec) 

(sec) 

(sec) 

1.00 

0.21 

0.17 

0.29 

0.25 

0.32 

0.28 

0.48 

0.1 

0.41 

0.44 

0.53 

0.075 

0.52 

0.51 

0.60 

0.0625 

0.59 

0.53 

0.63 

0.05 

0.82  ~ 

0.73  - 

0.69- 

0.025 

1.93  - 

>8- 

1.17- 

0.016 

>8- 

>8- 

2.66  - 

0.01 

>8- 

>8- 

>8- 

Figure  42.  Actuator  Failure  Identification  Performance  for  the  Standard 
Kalman  Filter  Bank  (SKFB)  with  a  Neyman-Pearson  Hypothesis  Testing 

Algorithm  (NPHTA). 
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4.4.4  Comparative  Failure  Identification  Performance.  We  used  the  failure 

identification  performance  data  shown  in  the  last  three  sections  to  compare  the  performance  for  the 
various  structures.  The  elevator  failure  identification  performance  of  all  three  structures  is  shown  in 
Figure  43.  Similarly,  the  performance  for  an  aileron  failure  is  shown  in  Figure  44,  and  the  rudder 
failure  identification  performance  is  shown  in  Figure  45. 

Note  that,  for  a  strong  input,  the  SKFB-NPHTA  structure  performs  the  best  of  the  three 
structures.  This  was  the  only  structure  that  did  not  show  a  decrease  in  failure  identification 
performance  for  stronger  inputs  for  the  rudder  failure.  This  structure  identified  the  failures  in  about 
half  the  time  that  it  took  for  the  SKFB-SHTA  structure,  and  about  one  tenth  the  time  it  took  for  the 
RCKFB-SHTA.  The  RCKFB-SHTA  structure  took  longer  at  these  input  levels  because  100  data 
samples  needed  to  be  collected  and  the  spectral  estimate  computed,  which  delayed  the  failure 
identification.  Collecting  a  small  number  of  data  samples  at  these  high  input  levels  might  produce 
faster  failure  identification  times. 

As  the  input  decreases,  the  failure  identification  performance  generally  decreases  for  the 
SKFB-SHTA  and  the  SKFB-NPHTA  structures,  and  seems  to  remain  fairly  constant  for  the 
RCKFB-SHTA  structure.  The  exception  to  this  general  rule  is  the  decrease  in  rudder  failure 
identification  performance  as  the  input  strength  increases  for  the  SKFB-SHTA  and  RCKFB-SHTA 
structures.  This  was  caused  by  both  the  relative  strength  of  the  rudder  input,  which  caused  a  larger 
residual  amplitude,  and  the  timing  of  the  failure,  which  occurred  just  as  the  residual  peaked,  so  the 
residual  took  longer  to  decay.  We  were  unable  to  conduct  a  failure  injection  timing  analysis  for  the 
RCKFB-SHTA  and  the  SKFB-NPHTA  structures,  but  we  believe  that  these  structures  would  suffer 
the  same  problems.  The  SKFB-NPHTA  does  not  show  this  behavior  on  these  plots  because  the 
failure  occurs  just  as  the  discrimination  measure  is  reset.  We  believe  that  the  SKFB-NPHTA  would 
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Figure  45.  Comparative  Failure  Identification  Performance  for  a 

Rudder  Failure. 


have  very  small  delays  due  to  the  timing  of  the  failure  injection.  These  conjectures  need  to  be 
researched  in  the  future. 

At  small  inputs  levels  the  performance  drops  of  rapidly,  especially  for  the  SKFB-NPHTA. 
The  RCKFB-SHTA  provides  the  best  failure  identification  performance  of  the  three  structures  at  loiv 
input  levels.  For  input  levels  below  1/20  of  the  original  input  used  for  this  research,  the  RCKFB- 
SHTA  structure  was  the  only  structure  that  still  provided  failure  identification  performance. 

All  three  structures  have  a  minimum  input  level  below  which  their  performance  diminishes 
rapidly.  As  the  input  strength  is  decreased,  the  approach  of  this  minimum  input  level  for  good  failure 
identification  performance  is  signaled  by  an  inaease  in  the  number  of  false  alarms.  This  suggests 
that  the  false  alarm  rate  could  be  monitored  (by  using  statistical  testing  for  estimating  the  false  alarm 
rate)  and  this  could  be  used  as  feedback  to  set  the  input  level  appropriately.  We  also  found  that  the 
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differences  in  failure  time  identification  for  the  SKFB-SHTA  could  be  minimized  by  setting  the  input 
level  just  above  the  point  where  the  failure  identification  performance  starts  to  deteriorate. 

These  results  strongly  suggest  that  a  combination  of  structures  would  provide  the  best 
failure  identification  performance,  This  combination  would  use  the  SKFB-NPHTA  structure  to 
provide  good  detection  at  strong  input  levels,  and  the  RCKFB-SHTA  structure  at  low  input  levels. 
Another  possibility  is  the  RCKFB  with  the  NPHTA,  which  may  provide  the  good  failure 
identification  performance  at  both  strong  and  low  input  levels.  Another  possible  structure  would 
combine  all  three  to  provide  redundancy,  which  could  then  be  used  to  decrease  the  false  alarm  rate  by 
using  voting  methods.  Also,  each  of  these  structures  could  be  implemented  using  equivalent 
residuals.  This  set  of  MMAE  structures  provides  a  rich  assortment  of  possible  combinations. 

4.5  Chapter  Summary 

In  this  chapter  we  have  presented  the  results  of  the  performance  evaluation  of  the  algorithm 
structures  that  were  developed  in  Chapter  3.  We  started  by  confirming  that  the  development  in 
Section  3.2,4  and  3.3.3,  where  we  characterized  the  mean  and  covariance  of  the  residual  in  the 
presence  of  mismodeling  in  the  input  control,  output,  and  state  transition  matrices.  We  corroborated 
this  characterization  for  mismodeling  in  input  control  and  output  matrices  by  making  graphical  and 
statistical  comparisons.  We  then  validated  the  development  of  the  equivalent  residuals  in  Section 
3.2.5.  We  compared  the  failure  detection  performance  for  three  different  MMAE  structures,  the 
Standard  Kalman  Filter  Bank  (SKFB)  with  the  Standard  Hypothesis  Testing  Algorithm  (SHTA),  the 
Residual  Correlation  Kalman  Filter  Bank  (RCKFB)  with  the  SHTA,  and  the  SKFB  (which  was  just  a 
"bank"  consisting  of  a  single  Kalman  filter)  with  the  Neyman-Pearson  Hypothesis  Testing  Algorithm 
(NPHTA).  We  found  that  the  SKFB-NPHTA  provided  the  quickest  failure  identification  at  fairly 


168 


strong  system  input  levels,  and  the  RCKFB-SHTA  provided  failure  identification  at  system  input 
levels  far  below  the  levels  required  by  the  other  two  structures. 
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y.  Conclusions 


5.1  Chapter  Overview 

In  this  chapter  we  present  our  conclusions  from  the  research  results  shown  in  Chapter  4.  In 
Section  5.2  we  look  at  the  characterization  of  the  mean  and  covariance  of  the  Kalman  filter  residual, 
in  the  presence  of  mismodeling  in  the  control  input  matrix,  output  matrix,  and  state  transition  matrix. 
We  present  our  conclusions,  in  Section  5.3,  concerning  an  equivalent  structure  to  the  Kalman  filter 
bank  in  which  the  residual  from  a  single  Kalman  filter,  along  with  its  state  estimates,  are  used  to 
compute  the  equivalent  residual  from  a  different  Kalman  filter.  In  Section  5.4,  we  draw  some 
conclusions  concerning  a  new  version  of  the  Kalman  filter  bank  in  which  spectral  estimates  of  the 
residuals  are  used  by  the  hypothesis  testing  algorithm  to  make  a  failure  status  evaluation.  Section 
5.5  presents  our  conclusions  for  the  Neyman-Pearson-based  hypothesis  testing  algorithm.  Finally, 
we  close  with  our  recommendations  for  future  research. 

5.2  Characterization  of  the  Residual 

We  developed  expressions  for  the  mean  and  covariance  of  the  Kalman  filter  residual  in  the 
presence  of  mismodeling  in  the  input  control  matrix,  output  matrix,  and  state  transition  matrix,  in 
Sections  3.2.3  and  3.2.4.  These  expressions  were  tested  by  comparing  the  computed  mean  and 
standard  deviation  to  actual  residuals  in  Section  4.2.  The  comparison  was  made  for  three  specific 
cases:  no  mismodeling,  a  mismodeled  control  input  matrix,  and  a  mismodeled  output  matrix.  We 
found  that  the  computed  mean  and  standard  deviation  accurately  characterized  the  actual  residual 
mean  and  standard  deviation,  thus  verifying  our  development. 
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5.3  Single  Residual  Kalman  Filter  Bank 

We  developed  a  methodology  in  Section  3.2.5,  in  which  the  residual  and  state  estimates  from 
one  Kalman  filter  can  be  used  to  compute  the  equivalent  residual  (and  state  estimates  if  needed)  from 
a  second  Kalman  filter,  if  the  system  model  used  by  the  second  filter  differs  from  the  first  filter's 
model,  only  in  the  input  control  matrix,  output  matrix,  and  state  transition  matrix.  We  propose  an 
MMAE-like  algorithm  composed  of  an  equivalent  Kalman  filter  bank  in  which  only  one  filter  is  fully 
implemented,  and  then  that  filter's  residual  and  state  estimate  are  used  to  compute  the  equivalent 
residuals  from  the  other  Kalman  filters  (that  are  not  actually  implemented)  in  the  filter  bank.  In 
Section  4.3,  we  compare  the  equivalent  residuals  with  the  actual  residuals,  and  found  that  they  are 
indeed,  equivalent.  This  structure  has  the  potential  to  attain  great  computational  savings,  since  only 
one  Kalman  filter  is  actually  implemented. 

We  noted  [39]  an  equivalence  class  relationship  of  this  structure  with  the  Generalized 
Likelihood  Ratio  (GLR)  failure  detection  structure.  The  GLR  failure  detection  structure  uses  a  set  of 
matched  filters  that  operate  on  the  residual  from  a  single  Kalman  filter.  Each  of  the  matched  filters  is 
based  on  a  hypothesized  behavior  of  the  residual,  such  as  a  sudden  additive  bias,  a  slowly  increasing 
bias,  an  increase  in  measurement  noise  covariance,  or  no  change  in  the  residual.  The  hypothesized 
residual  behaviors  are  based  on  the  designer's  expectations  that  a  certain  failure  will  cause  this 
behavior  in  the  residual,  and  they  are  often  rather  ad  hoc  and  simple  in  practice.  This  interpretation 
of  our  results  holds  that  the  linear  transforms  operating  on  the  output  of  a  single  Kalman  filter  (as 
seen  in  Figure  6, 25,  and  27)  are  equivalent  to  the  best  possible  "matched  filters"  of  the  GLR 
structure.  This  allows  the  designer  to  postulate  a  certain  failure  model  and  compute  the  linear 
transform  for  that  model  using  the  postulated  differences  between  the  fully  functional  model  and  the 
hypothesized  failure  model.  The  linear  transform  is  equivalent  to  the  best  possible  GLR  matched 
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filter  for  that  failure  model.  The  interpretation  establishes  an  equivalence  class  relationship  between 
MMAE  failme  identification  and  the  hest-possihle  GLR  failure  detection  and  isolation  systems. 

5.4  Residual  Correlation  Kalman  Filter  Bank 

We  developed  a  new  kind  of  filter  bank  structure  in  Section  3.2.6,  in  which  an  estimation  of 
the  spectral  content  of  the  residuals  (or  equivalent  residuals)  is  used  by  the  hypothesis  testing 
algorithm.  The  actuator  failure  identification  performance  of  this  Residual  Correlation  Kalman  Filter 
Bank  (RCKFB)  was  tested  and  the  results  were  presented  in  Section  4.4.2,  and  a  comparison  of  this 
performance  against  other  techniques  was  made  in  Section  4.4.4.  In  Section  3.2.4,  we  showed  that, 
when  an  actuator  failure  occurs,  elements  of  the  control  input  are  added  to  the  residual.  If  the  control 
input  elements  are  sinusoids,  then  the  residual  includes  a  sinusoid,  with  the  same  frequency  as  the 
control  input  elements.  Therefore,  we  can  use  the  spectral  content  of  the  residual,  at  the  input 
frequencies,  to  indicated  the  presence  of  an  actuator  failure.  We  found  that,  when  the  control  input  is 
fairly  strong,  this  technique  takes  longer  than  the  standard  MMAE  to  make  the  correct  failure 
identification,  due  to  the  time  used  to  collect  a  sufficient  number  of  samples  to  generate  an  adequate 
spectral  estimate.  However,  this  technique  provides  failure  identification  performance  at  small  input 
levels  where  the  standard  MMAE  no  longer  functions.  This  technique  looks  quite  promising  as  a 
method  to  provide  good  failure  identification  with  subliminal  input  dithers,  which  could  be  used  in 
conjunction  with  the  standard  MMAE  structure  to  provide  failure  identification  performance  at  low 
input  levels  and  to  corroborate  the  failure  identification  of  the  SKFB-SHTA  structure  at  higher  input 
levels. 
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5.5  Nevman-Pearson  Hypothesis  Testing  Algorithm 


We  developed  a  Neyman-Pearson  Hypothesis  Testing  Algorithm  (NPHTA)  in  Section  3.3.2. 
We  first  extended  the  standard  development  for  a  scalar  process,  to  a  multidimensional  process,  and 
then  extended  that  from  a  single  time  sample  hypothesis  test,  to  multiple  time  sample  hypothesis 
testing.  Sequential  Probability  Ratio  Tests  (SPRT)  [34]  have  been  shown  to  be  superior  to  the 
Neyman-Pearson  most  powerful  tests  when  operating  on  sequential  data,  because  the  SPRT  allows 
variable  sample  size  and  makes  a  decision  when  and  only  when  sufficient  information  is  obtained. 
We  extended  the  NPHTA  to  multiple  time  samples  by  incorporating  a  triggering  methodology,  and 
so  our  algorithm  should  share  many  of  the  attributes  of  the  SPRT  algorithm. 

We  showed  that,  by  using  the  NPHTA,  we  could  perform  multiple  hypothesis  tests  on  a 
single  Kalman  filter  residual,  instead  of  having  to  form  a  residual  for  each  hypothesis.  Thus,  we 
would  use  a  Kalman  filter  "bank"  of  a  single  Kalman  filter,  followed  by  the  NPHTA.  This  structure 
was  implemented  and  tested,  and  the  results  were  presented  in  Section  4.4.3.  The  actuator  failure 
identification  performance  was  compared  to  other  MMAE  structures  in  Section  4.4.4.  We  found  that 
this  technique  was  significantly  faster  than  the  standard  MMAE,  when  the  input  strength  was  large 
enough  to  provide  good  failure  identification  performance.  When  the  input  strength  dropped  below 
this  minimum  level  that  was  needed  for  good  failure  identification,  the  performance  of  this  technique 
dropped  off  dramatically.  The  performance  for  sensor  failures  was  found  to  be  quite  good  for 
sensors  that  do  not  have  strong  measurement  noise,  but  the  false  alarm  rate  was  rather  high  for 
sensors  that  have  strong  measurement  noise.  The  performance  of  the  NPHTA,  along  with  the  single 
Kalman  filter  "bank"  structure,  argue  strongly  for  implementing  this  structure  for  applications  for 
which  the  input  strength  can  be  above  the  minimum  level  needed  for  good  failure  identification. 
However,  if  the  input  levels  are  below  this  minimum  level,  the  RCKFB  structure  would  need  to  be 
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implemented.  Thus,  we  foresee  having  a  two-structure  MMAE  algorithm;  a  RCKFB-based  MMAE 
for  low  input  strengths,  and  an  NPHTA  for  high  input  strengths.  A  three-structure  MMAE  may  be 
desired,  where  the  third  structure  is  the  standard  MMAE,  to  corroborate  the  results  from  the  other 
two  structures.  All  of  these  structures  could  use  the  equivalent  Kalman  filter  bank  composed  of  a 
single  Kalman  filter  with  a  series  of  linear  transforms  to  produce  the  equivalent  residual  of  a  fully 
implemented  SKFB. 

5.6  Recommendation.s  for  Future  Research 

This  investigation  opens  up  several  possible  avenues  of  further  research.  First,  and 
foremost,  would  be  combining  the  RCKFB  with  the  NPHTA.  This  has  the  potential  to  provide  the 
excellent  failure  identification  performance  of  the  NPHTA  at  high  input  strength,  with  the  good 
failure  identification  performance  of  the  RCKFB  at  low  input  strength.  Along  with  this,  the 
performance  of  the  RCKFB  for  various  values  of  N,  the  number  of  data  points  used  by  the  spectral 
estimator,  needs  to  characterized.  Further,  the  use  of  other  methods  of  spectral  estimation,  such  as 
the  Blackman-Tukey,  model-based  methods,  and  Maximum  Entropy  Methods  (MEM),  needs  to  be 
investigated.  The  failure  identification  performance  of  these  other  techniques  should  be  evaluated  for 
each  of  the  different  methods,  and  a  comparison  made  with  the  periodogram  that  was  used  for  this 
research.  Also,  the  performance  of  the  NPHTA  needs  to  be  explored  for  various  detection  and  false 
alarm  probability  settings,  and  a  method  of  determining  the  NPHTA  sensor  failure  performance 
needs  to  be  developed. 

Several  of  the  MMAE  structures  that  were  developed  in  the  research,  show  great  promise  in 
reducing  the  computational  loading  of  the  MMAE:  specifically,  the  single  residual  Kalman  filter 
bank,  developed  in  Section  3.2.5,  and  the  single  Kalman  filter  bank  used  by  the  NPHTA.  The 
computational  savings  attained  by  these  structures  need  to  be  determined,  for  various  applications. 
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by  further  research.  It  is  particularly  attractive  for  actuator  failure  identification  because  of  the 
sparseness  of  the  system  input  matrix  that  is  used  by  the  linear  transformations  that  replace  the  many 
Kalman  filters  of  a  standard  MMAE  algorithm.  Implementing  this  structure  for  sensor  failure 
identification  may  produce  computational  savings  if  the  difference  between  the  Kalman  filter  gain 
matrices  results  in  a  sparse  matrix.  Also,  the  computational  cost  of  implementing  the  RCKFB  using 
Eq  (71),  in  which  the  full  estimation  of  the  autocorrelation  matrix  is  used,  needs  to  be  compared  to 
the  savings  attained  by  implementing  Eq  (72),  where  only  the  diagonal  terms  of  the  estimated 
autocorrelation  matrix  are  used,  thus  ignoring  the  aoss-correlation  terms.  This  cost  savings  needs  to 
be  determined  and  compared  to  any  reduction  in  failure  identification  performance  caused  by 
ignoring  the  cross-correlation  terms. 

The  variation  in  failure  identification  times  caused  by  different  failure  injection  times  needs 
to  be  examined.  The  design  of  "balanced"  actuator  inputs,  in  which  the  input  for  one  actuator  has  an 
equivalent  effect  on  the  likelihood  quotient  as  the  input  for  another  actuator,  needs  to  be  a  part  of  this 
examination.  This  would  help  determine  the  appropriate  mix  of  inputs  to  provide  good  failure 
identification  of  aU  pertinent  failure  modes. 

This  suggests  an  alternative  structure  to  the  RCKFB,  one  that  uses  the  spectral  content  of  the 
scalar  likelihood  quotients  from  the  various  Kalman  filters,  rather  than  the  spectral  content  of  the 
residual  vector.  This  structure  would  exploit  the  cross-coupling  terms  that  were  discarded  when  we 
implemented  the  Fourier  transform-based  spectral  estimator.  Thus,  this  structure  would  use  the 
simpler  and  less  computationally  intensive  Fourier  transform-based  estimator,  as  in  Eq  (72),  but  still 
include  all  of  the  cross-terms  that  are  in  the  periodogram-based  spectral  estimator  of  Eq  (71).  The 
development  of  the  algorithm  shows  great  promise  for  improving  the  RCKFB  failure  identification 
performance. 
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Future  research  should  aim  at  an  information-theoretic  development  of  a  class  of  algorithms 
that  fully  (or  partially)  exploits  the  information  content  in  the  residuals.  The  algorithm  structures 
that  were  developed  in  this  research  would  be  examples  of  elements  of  such  a  class  of  algorithms. 

A  fuller,  more  rigorous  investigation  of  the  interrelationship  between  the  equivalent-residual 
MMAE-like  structures  and  the  GLR  structure  that  uses  the  best-possible  matched  filters  merits  being 
accomplished.  The  linear  transforms  of  this  research  are  the  best  matched  filters  in  the  sense  that 
they  truly  represent  the  time  history  of  the  residuals  caused  by  a  physically  motivated  failure  model. 
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